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Abstract 

In this work we investigate the relationship among three classical sampUng techniques: the inverse 
of density (Khintchine's theorem), the transformed rejection (TR) and the generaUzed ratio of uniforms 
(GRoU). Given a monotonic probabiUty density function (PDF), we show that the transformed area 
obtained using the generaUzed ratio of uniforms method can be found equivalently by applying the 
transformed rejection sampUng approach to the inverse function of the target density. Then we provide 
an extension of the classical inverse of density idea, showing that it is completely equivalent to the GRoU 
method for monotonic densities. Although we concentrate on monotonic probability density fimctions 
(PDFs), we also discuss how the results presented here can be extended to any non-monotonic PDF that 
can be decomposed into a collection of intervals where it is monotonically increasing or decreasing. In 
this general case, we show the connections with transformations of certain random variables and the 
generalized inverse PDF with the GRoU technique. Finally, we also introduce a GRoU technique to 
handle unbounded target densities. 

Index Terms 

Transformed rejection sampling; inverse of density method; Khintchine's theorem; generalized ratio 
of uniforms technique; vertical density representation. 
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I. Introduction 

Monte Carlo (MC) methods are often used for the implementation of optimal Bayesian estimators in 



many practical applications, ranging from statistical physics |Rosenbluth and Rosenbluth 19551 Siepmann 



and Frenkel[ 1992 1 to nuclear medicine [Ljungberg et al.[ 1998 1 and statistical signal processing |Djuric 



et aH |2003[ [Martino and Mfguez] |2010[ [Ruanaidh and Fitzgerald[ |1996[ . Many Monte Carlo techniques 



have been proposed for solving this kind of problems either sequentially (SMC methods, also known as 
particle filters), making use of Markov chains (MCMC methods) or otherwise [Fitzgerald 2001 Gilks 



et aH 1 1995] [Ruanaidh and Fitzgerald[|1996| |. Sampling techniques (see e.g. |Devroye[ |1986[ [Gentle[ |2004 



Hormann et al. 2003 1 for a review) are the core of Monte Carlo simulations, since all of them rely on the 



efficient generation of samples from some proposal PDF |Liu 2004[ Robert and Casella) [2004 1. Many 
sampling algorithms have been proposed, but the problem of drawing samples efficiently from a generic 
distribution is far from trivial and many open questions still remain. 

In this paper we investigate the relationship between three classical sampUng techniques: 

• the inverse-of-density technique | Devroye[ 1986 Khintchine 1938 1, 



• the transformed rejection sampling | |von Neumann} |1951[ |Wallace[ |1976[ , also known as "almost 
exact inversion" | Devroye| 1986[ Chapter 3] and "exact approximation method" [MarsagUa 1984 1, 

• and the ratio-of-uniforms method [Kinderman and Monahan 1977[ Wakefield et al.[ 1991 1. 

We present new connections among them, useful to design more efficient sampling techniques. Although 
in the sequel we concentrate mainly on monotonic PDFs, we also discuss the relationships among these 
techniques in more general cases, especially in the last two sections and the Appendix. 

The first method considered, the inverse-of-density (loD) technique | Devroyej 1986 Chapter 4], 
IDevroye 1984 Isii 1958 Jones 2002 Khintchinel 1938[ (often known as Khintchine' s theorem [Feller 



T97T] pp. 157-159], pryson and Johnsonj [T9821 [Chaubey et alj [20T0l [Jones] [20021 [Khintchine[ [T938 



Olshen and Savage] [1970[ [Shepp[ |1962[ , both for monotonic PDFs and for symmetric unimodal PDFs), 



is a classical sampling technique. Given a monotonic target PDF, pq{x) = Kp{x) (also denoted often 
as po{x) oc p{x), omitting the normalization constant, K > 0), this method provides a closed-form 
relationship between the samples from the PDF defined by the unnormalized inverse density, p^^{y), 
and the desired samples, distributed according to the normalized PDF, pq{x). Hence, if we are able to 
draw samples easily from p^^iy), then it is straightforward to generate samples from the target PDF by 
using the loD approach. Clearly, the practical applicability of the loD method depends on the feasibility 
of drawing samples from the inverse density p~^{y). 
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The loD method can be easily extended to non-monotonic densities (see e.g. |Devroye 1986 Jones 



2002| |) both unidimensional and multidimensional | |Bryson and Johnson 1982[ de Silva[ 1978[ . Moreover, 



the loD presents several relationships (see e.g. | |Jones[[2002| ) with vertical density representation (VDR) 
I Fang et al.| [200T] |Kotz et al.| [19971 |Kotz and Trouttl [T996| |Kozubowski| [20021 [Troutt[ \m\\ [T993 



Troutt et"aLl [2004[ , especially with the so-called second type VDR |Fang et all [200 , |Troutt et~ar 



20041 Chapter 3]. 

The second tackled method is transformed rejection sampling (TRS) | Devroye[|1986[ Marsaglia 1984 



Wallace! 1976[ . The rejection sampling (RS) is another standard Monte Carlo technique that use a a 



simpler proposal distribution,7r(x) to generate samples and, then, to accept or discard them according to 
a ratio between the target and proposal densities -^^^ (where Lir{x) > p{x)). Hence, the fundamental 
figure of merit of a rejection sampler is the mean acceptance rate (i.e. the expected number of accepted 
samples out of the total number of proposed candidates). 

The most favorable scenario for using the RS algorithm occurs when the unnormalized target PDF, 
p{x), is bounded with bounded domain. In this case, the proposal PDF 7r(x) can be a uniform density (the 
easiest possible proposal), and calculating the bound L for the ratio p(x)/7r(x) is equivalent to finding an 



upper bound for the unnormalized target PDF, p{x), which is in general a much easier task | Devroye 1986 



Hermann et al. 2003 1. Indeed, in this scenario several sophisticated and efficient acceptance/rejection 



methods that achieve a high acceptance rate have been devised: adaptive schemes | Gilks and Wild[ 1992 



Martino and Mfguez 2011a|, strips techniques [Hormann et al.[ 2003 [ Chapter 5], | Ahrens 1993 1995 



Marsaglia and Tsangl [2000| , patchwork algorithms [Kemp[ [1990[ [Stadlober and Zechner} |1999[ |, etc. 



However, in general the target p{x) can be unbounded or with an infinite support and the choice 



of a good proposal PDF becomes more critical (see, for instance Martino and Mfguez |2011b|). 
In order to overcome this problem, different methods have been proposed to transform the region 
corresponding to the area below p{x) into an alternative bounded region. A straightforward solution from 



a theoretical point of view is the transformed rejection sampling (TRS) [Botts et al. 2011 Hormann 



1993[ Hormann and Derflinger[ [1994[ [Marsaglial 1984[ Wallace 1976|, which is based on finding a 



suitable invertible transformation, f{x) : Vx — T^z, such that the region below p{x) is transformed into 
an appropriate bounded set. Making use of this transformation we can define a random variable (RV) 
Z = f{X), with unnormalized PDF p{z) = p{f~^{z)) x and denoting the derivative of 

f^^{z), draw samples {z^^\ . . . ,z^^^} from p{z), and convert them into samples from the target PDF, 
{x^^\ . . . , x^^)} = {f^^{z^^^), • • • , f^^{z^^^)} by inverting the transformation. 

Obviously, attaining a bounded PDF p{z) requires imposing some restrictions on the transformation 
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f{x) that depend on the unnormaUzed target PDF, p{x) [Hormann and Derflinger 1994 Wallace 



|1976| |. Furthermore, the performance of the TRS approach depends critically on a suitable choice of 
the transformation function, f{x). Indeed, if f{x) is chosen to be similar to the cumulative distribution 
function (CDF) of the target RV, Fx{x), the PDF of the transformed RV, p{z), becomes flatter and 
closer to a uniform PDF and higher acceptance rates can be achieved. In particular, if f{x) = Fx{x), 
then p{z) is the uniform density in [0,1], implying that we can easily draw samples from it without 
any rejection and justifying the fact that this technique is sometimes also called almost exact inversion 
method | |Devroye[ |1986| . 

Another approach to work with bounded region is the so-called ratio-of-uniforms (RoU) technique 
I Devroye| 1986[ Kinderman and Monahan 1977[ (the third technique that we address here). The RoU 
ensures that, given a pair or independent RVs, {V,U), uniformly distributed inside Ar = {{v,u) G : 
< u < \/p{v/u)}, then x = v/u i?, distributed exactly according to the target PDF, po{x). Hence, in 
the cases of interest (i.e. when the region Ar is bounded) the RoU provides us with a bidimensional 
region, Ar, such that drawing samples from the univariate target density is equivalent to drawing samples 



uniformly inside Ar, which can be done efficiently by means of rejection sampling schemes |Luengo 



and Martinol |2012[ |Leydold[ |2000[ |2003[ [Perez et al.] |2008| . Unfortunately the region Ar provided by 
RoU is only bounded when the tails of the target density decay faster than which is not always 

fulfilled for the PDFs of interest. 

Consequently, several generalizations of the RoU method have been proposed in the literature (see 



e.g. IJones and Lunn 1996 Wakefield et al. 1991 1 and more related materials that can be found in 
Parbul [T9821 |Curtiss| [T94T] |Dieter[ [T9891 piiaglial [T965| [Perez et al.| [2008| [Stefanescu and Vaduva 



1987[ [Vaduva[ [1982[ ). The most popular of those extensions is the so called generalized ratio-of- 
uniforms (GRoU) [Wakefield et al. 1991 1, which shows that x = v/g{u) is distributed according 
to the target PDF, pq{x), when the random vector {V, U) is uniformly distributed inside the region 
Ag = {{v,u) G : < u < g^^{cp{v / g{u)))}, with c > being a constant term and g{u) a strictly 
increasing differentiable function on such that g{{)) = 0. 

These two techniques (TRS and GRoU) have been introduced separately in the literature and their 
connection has not been explored as far as we know. The primary goal of this paper is showing that 
there is a close relationship between both approaches. Indeed, one of the main results in this work is 
proving that the transformed region attained using the GRoU technique [ [Wakefield et al. 1991 1 can also 
be obtained applying the transformed rejection approach |Wallace[ 1976 1 to the unnormalized inverse 
PDF, p^^{y), for monotonic target PDFs, po{x) oc p{x). Moreover, we introduce an extended version 
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of the standard inverse-of-density method |Devroye 1986[ Chapter 4], | Jones[ 2002[ [Khintchine[ 1938 1, 



which is strictly related to the GRoU method and show that the GRoU sampling technique coincides 
with this extended version of the inverse-of-density method. 

Considering a monotonic unnormalized target PDF, p{x), in this work we show that the region Ag 
defined by the GRoU can be obtained transforming an RV Y with unnormalized inverse PDF p^^{y), 
and that the relationship between the points in this region Ag and the samples drawn from p{x) is 
provided by the novel extended version of the loD method, introduced here. Hence, as a conclusion we 
can assert that, for monotonic PDFs, the GRoU can be seen as a combination of the transformed rejection 
sampling method applied to the unnormalized inverse PDF, p^^(y), and an extended inverse-of-density 
technique. We also investigate the connections among TRS, loD and GRoU for generic non-monotonic 
target PDFs. Finally, taking advantage of the previous considerations we introduce a GRoU technique to 
handle unbounded target distributions. 

The rest of the paper is organized as follows. In Section [II] we provide some important considerations 
about the notation, we formulate the fundamental theorem of simulation (which is the basis for all the 
sampling methods discussed), and briefly describe the standard inverse-of-density and rejection sampling 



techniques, thus providing the background for the rest of the paper. Then, Sections [III] and |IV] provide 
a detailed description of the two sampling methods compared, transformed rejection sampling and the 
generalized ratio-of-uniforms respectively, focusing on the different possible situations that may be found 
and particularly on the conditions required for obtaining finite sampling regions. This is followed by 



Section W\ where we introduce an extension of the inverse-of-density method, and Section VI which 



provides the main result of the paper: the relationship between the ratio-of-uniforms, transformed rejection 



and the inverse-of-density methods. Section [Vll provides some further considerations about the different 



approaches considered, whereas Section VIII discusses their extension to non-monotonic PDFs. Section 



IX is devoted to design a GRoU for unbounded distributions, using the previous considerations and 



observations. Finally, the conclusions and the appendix close the paper. 

II. Background 

A. Important consideration about the notation 

In the sequel we always work with proper but unnormalized PDFs, meaning that integrating them over 
their whole domain results in a finite positive constant, but not necessarily equal to one. As an example, 
consider the normalized target PDF, pq{x) = Kp{x), with K > Q denoting the normalization constant. 
All the subsequent methods will be formulated in terms of p{x), which is the unnormalized target PDF, 
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Since 

1 



p{x)dx = — , (1) 

with K > 0, but K I in general. Hence, the integral is finite but not necessarily equal to one. 

Furthermore, in order to get rid of the normalization constant, we will also work with the unnormaUzed 
inverse target PDF, p^^{y), for monotonic target PDFs or its generalized version, p^{y), for non- 
monotonic PDFs. Note that the normalized inverse target PDF, pQ^{y), can no longer be obtained from 
the unnormalized inverse target PDF, p^^{y), simply multiplying by a normalization constant. A scaling 
of the independent variable, y, must be performed instead in order to attain PQ^{y) = p^^{y/K). We 
remark also that Kp^^{y) ^ Po^(y) = P^^iu/K), i.e. the normalized version of the unnormalized 
inverse target PDF will be different from the normalized inverse target PDF in general. This is due to the 
fact that the support of ^(y) will usually be different from the support of p^^{y), due to the scaling 
of the independent variable, y, performed on p^^(y) in order to obtain Pq ^(y). Finally, note also that, 
given a sample y' from Kp~^{y) we can easily obtain samples from the normalized inverse target RV, 
PQ^{y), since y' /K PQ^{y)- AH these issues ai^e clearly illustrated in the following example. 

Example 1: Consider a half Gaussian random variable with the following PDF: 



Po(:r) = y^exp(-^), (2) 

for < a; < oo. The half Gaussian PDF given by Q is bounded, < po{x) < ^J2f{j^a^, with bounded 
support, T>x = = [0, oo), and we can easily identify the corresponding unnormalized PDF, 

p{x) = exp (^-^^ ' (3) 

for < X < oo, and the normalization constant, 

K = \\ (4) 

The functional inverse of (|2]) (i.e. the normalized inverse taget PDF) is given by 



Po^(2/) = A/^^log — (5) 

for < y < \fljjjuj^, and log indicating the natural logarithm. Note that Q defines a proper normalized 
PDF, since Pg ^(y) > for any value of y and 

p-^\y)Ay = \. (6) 
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Furthermore, it is an unbounded PDF, since limj^^j^o Pq {u) = oo, but with a bounded support, 
Vy = (0, ^y2/{7^a'^)]. Similarly, the unnormalized inverse target PDF is given by 

p-\y) = V-2a2log y, (7) 

for < y < 1, which is also a proper unbounded PDF with a bounded support, Vy = (0, 1]. We note 
that the normalized version of (|7]l, Kp^^{y) = 2^ —^\og y for < y < 1, is clearly different from the 
normalized inverse target PDF, ^(y) = p^^{y/K) for < y < ^J^jijxa'^), given by Q. Finally, we 
also notice that, given y' ~ Kp~^{y), then y' /K = y' ^ -na'^ /2 is distributed as p^^(y/K) = PQ^{y), as 
discussed before. 

In order to conclude this section, it is important to remark that all the discussions and algorithms 
shown below do not require the knowledge of the normalization constant. Hence, we can work with 
unnormalized PDFs without any loss of generality, since all the results attained in the sequel can also be 
formulated using normalized PDFs, although the notation becomes more cumbersome. This is a standard 
approach, followed not only by the RS and TRS methods, but also by most other standard sampling 
algorithms, like the RoU or the GROU. Therefore, in the sequel we will use X ~ p{x) and Y ~ p^^{y) 
to indicate that the PDFs of the RVs X and Y are proportional to the unnormalized target and inverse 
target RVs, even though p{x) and p^^{y) will not be normalized in general. See the Appendix for a 
detailed revision of the notation used throughout the paper 

B. Fundamental theorem of simulation 

Many Monte Carlo techniques (inverse of density method, rejection sampling, slice sampling, etc.) 
are based on a simple result, known as the fundamental theorem of simulation, that we enunciate in the 
sequel. 

Theorem 1: Drawing samples from a unidimensional random variable X with probability density 
function pq{x) = Kp{x), where K > is a constant, is equivalent to sampling uniformly inside the 
bidimensional region 

A = {(x,y)GM': Q<y<p{x)}. (8) 



Proof 1: Straightforward. See [ [Robert and Casella| [20041 Chapter 2]. 

Hence, according to Theorem [TJ if the pair of random variables (X, Y) is uniformly distributed inside 
the region Aq, which corresponds to the area below p{x), then the PDF of X is proportional to p{x). 
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whereas the random variable Y plays the role of an auxiliary variable. Many Monte Carlo techniques 
make use of this theorem explicitly to simulate jointly the random variables (X,Y), discarding Y 
and considering only X, which is a univariate random variable marginally distributed according to 
the unnormalized target PDF, p{x) [Robert and Casella 2004 1. Figure [T] depicts an example of an 
unnormalized target PDF, p{x), and the region delimited by it.The two methods described in the 
sequel, inverse of density and rejection sampling, are clear examples of how this simple idea can be 
applied in practice to design Monte Carlo sampling algorithms. 




Fig. 1. The region corresponding to the area below the unnormalized target PDF, p{x). 



C. Inverse of density method for monotone PDFs 
In this section we present the inverse of density (loD) method [Devroye 1986[ Chapter 4], | Jones 



20021, often known as Khintchine's theorem | Feller[ 1971 pp. 157-159], [Bryson and Johnson 1982 



de Silval [T978| [Ml [T9581 |Khintchine| [T9381 [Olshen and Savage| [T970| , both for monotonic and for 



symmetric unimodal densities | Chaubey et al. 2010[ Shepp[ [1962 1 . Note once more that, although we 
concentrate here on monotonic PDFs, this result can be easily extended to generic PDFs, as tackled in 



Section VIII- A and also shown in | Devroye 1986 1. The standard formulation for the loD method is 
the following. Let us consider a monotonic unnormalized target PDF, p{x), and denote by p^^{y) the 
corresponding inverse function of the unnormalized target density. 

The fundamental idea underlying the loD approach is noticing that p~^{y) can also be used to describe 
^o> as illustrated graphically in Figure [2] Consequently, the region associated to p{x), 

A = {(x,2/)eM': 0<y<p{x)}, (9) 

shown in Figure 2(a)[ can be expressed alternatively in terms of the inverse PDF as 

A = {(y,x)GM2. o<x<p-i(y)}, (10) 
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as depicted in Figure 2(b) Therefore, we can proceed in two alternative ways in order to generate samples 
{x',y') uniformly distributed inside Aq- 
1) Draw first x' from p{x) and then y' uniformly in the interval [0,^(2;')], i.e. y' ~ U{[0,p{x')]), as 



shown in Figure 2(a) 



2) Draw first y' from p ^{y) and then x' uniformly in the interval [0,p ^{y')], i.e. x' L({[0,p ^{y')]), 
as shown in Figure |2(b)[ 

Both procedures allow us to generate points (x', y') uniformly distributed inside the region Aq- Moreover, 
from the fundamental theorem of simulation, the first coordinate x' is distributed according to the 
unnormalized target PDF, p{x), whereas the PDF of the second coordinate y' is proportional to the 
unnormalized inverse PDF, p^^{y). Hence, the key idea of the inverse of density method is that, whenever 
we are able to draw samples y' from p^^{y) more easily than samples x' from p{x), we can use the 
second procedure to generate samples x' from p{x) more efficiently. 




(a) (b) 

Fig. 2. Two ways of drawing a random point {x',y') uniformly inside the region Ao- (a) Draw first x' from p{x) and then 
y' ~ W([0,p(s')]). (b) Draw first y' from p'^{y) and then a;' ~ W([0,p"^(i/')])- 



Note that generating a sample x' uniformly inside the interval [0,a], i.e. x' ~ U{[0,a\), is equivalent 
to drawing a sample w' uniformly inside [0, 1] and then multiplying it by a, i.e. x' = w'a. Thus, given a 
known value y', drawing a sample x' uniformly inside the interval [0,p^^{y')], i.e. x' ~ U{[0,p^^{y')]), 

'Noting that the samples y' generated in this way are distributed according to p~^{y), we remark that this method can 
always be used to generate samples from the generalized unnormalized inverse PDF, pQ^{y), even when p{x) is non-monotonic. 
However, in this case the geometric interpretation of this generalized inverse PDF becomes more complicated, since its definition 



may not be straightforward, as shown in Section VIII 
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is equivalent to generating a sample w' uniformly inside [0, 1] and then taking 

/ / -1/ /\ 



(11) 



which is the expression frequently provided for the loD method. We also remark that, for a proper 
monotonia unnormalized density p{x), p~^{x) is also a proper monotonia unnormalized PDF, obtained 
simply through functional inversion. 

Obviously, the interest in using this technique depends on the feasibility of drawing samples from the 
unnormalized inverse PDF, p^^{y), more easily than from the unnormalized target PDF, p{x), as already 
mentioned. The following example shows a practical application where the loD method provides a clear 
advantage over the direct generation of a random variable. 

Example 2: Assume that we need to draw samples from 



Po{x) = p{x) 



'log 



TTX^ 



with < a; < 



(12) 



Since p^^{y) = Y|exp(— y^) = exp(— y^/2) is the half Gaussian PDF used in the previous 

example, we can easily draw y' from p^^{y), then w' from a uniform PDF inside [0, 1], and finally 
obtain a sample x' = w'p^^{y'), which is distributed according to the target PDF, p{x). 

Finally, we notice that it is possible to find this method in other forms related to vertical density 
representation in the literature [Fang et aL| |2001[ |Jones[ |2002[ IKhintchine] |1938[ [Troutt et aLj |2004| |. 
Indeed, let us consider a random variable Y which follows a strictly decreasing unnormalized PDF, 
p^^{y)- Then, the random variable U = p~'^{Y) is distributed as 



l{u) = -p {p{u)) 



dp{u) 



_ dpiu) 
-u- 



(13) 



du du 

and this unnormalized PDF, q{u), is called the vertical density w.r.t. p~^{y). Making use of this result, 
the inverse of density method, summarized by equation ( [TT| ), can be expressed alternatively in this way: 
given w' ~ ^([0, 1]) and u! ~ q{u), then the sample 

x' = w'u! = w'p-^{y'), (14) 



is distributed as p{x) provided that y' is a sample from p ^(y). The relationship in Eq. (14i is usually 
known as Khintchine's theorem. 



D. Rejection sampling 

Another technique that clearly applies the simple idea exposed in Section II-B is rejection sampling. 
Rejection sampling (RS) is a universal method for drawing independent samples from an unnormalized 
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target density, p{x), known up to a proportionality constant K > 0. Let ■k(x) be a (possibly unnormalized) 
proposal PDF and L an upper bound for the ratio p{x)/tt{x), i.e. 

L>'-^. (15) 

RS works by generating samples from the proposal PDF, 7r{x), and accepting or rejecting them on the 
basis of this ratio. The standard RS algorithm can be outlined as follows. 

1) Draw x' ~ tt{x) and w' ~ ^^([0, 1]). 

2) If w' < j^^, then x' is accepted. Otherwise, x' is discarded. 

3) Repeat steps 1-2 until as many samples as required have been obtained from the target PDF. 
Alternatively, the procedure undertaken by the RS method can also be summarized in the following 
equivalent way that remarks its close connection to the fundamental theorem of simulation. 

1) Draw x' from tt{x). 

2) Generate y' uniformly inside the interval [0,L7r(x')], i.e. y' ~ ZY([0, L7r(x')]). 

3) If the point {x',y') belongs to ^o> the region corresponding to the area below the unnormalized 
target PDF p{x) as defined by ([8]l, the sample x' is accepted. 

4) Otherwise, i.e. whenever the point (x', y') falls inside the region located between the functions L7r{x) 
and p{x), the sample x' is rejected. 

5) Repeat steps 1-4 until as many samples as required have been obtained from the target PDF. 
Figure |3] provides a graphical representation of the rejection sampling technique. Here, the green region 
corresponds to as defined by ([8]l, the region associated to the target PDF inside which we want 
to sample uniformly (i.e. the acceptance region), whereas the red region indicates the region located 
between the functions L'k{x) and p{x), where we do not want our samples to lie (i.e. the rejection 
region). Defining 

X = {{x, y) e : < y < (16) 

this rejection or exclusion region is given by the set-theoretic difference or relative complement of 
inside A-j^: 

AI=A^\Aq = {{x, y) G M2 . ^ y < L7r(x)}. (17) 

Now, the RS algorithm proceeds by drawing first a sample from the proposal PDF, x' ~ vr(x), and then a 
second sample from a uniform distribution, y' ~ ^([0, Lt\:{x')]). If the point (x', y') belongs to (green 
region), as it happens for the point indicated by a filled dark green circle in Figure [5] the sample x' is 
accepted. Otherwise, whenever the point {x\y') belongs to Aq (red region), as it happens for the point 



12 



indicated by a filled dark red circle in Figure |3j it is discarded. Note that, since y' ~ ^([0, L7r(x')]) can 
be expressed alternatively as y' = Ltt{x')w' with w' ~ ^([0, 1]), the previous conditions are equivalent 
to accepting x' whenever y' = Ltt{x')w' < p{x'), which happens if and only if {x',y') belongs to the 
green region, and rejecting x' otherwise, i.e. whenever y' = Ltt{x')w' > p{x'), which happens if and 
only if {x',y') belongs to the red region. This is equivalent to the condition shown in step 2 of the first 
formulation, demonstrating the equivalence between both descriptions of the RS algorithm. 




Fig. 3. Graphical description of the RS procedure. The green region corresponds to the acceptance region Ao, as defined by 
([8|, whereas the red region indicates the rejection region, Ao, located between the functions Ln{x) and p{x), as defined by 
^IT^, and x' denotes a sample drawn from the proposal PDF, n{x), for the two possible situations that can occur: x' £ 
(filled dark green circle) and x' £ Aq (filled dark red circle). 



The fundamental figure of merit of a rejection sampler is the mean acceptance rate, i.e. the expected 
number of accepted samples out of the total number of proposed candidates, which is given by 



I'^vrl I'^ol ~l~ 1-^0 1 |-^7r| 



where \C\ denotes the Lebesgue measure of set C, and the last two expressions arise from the fact that, 
since A-k = ^qU^q ^i^^ ^qH^q = 0, then |^7r| = |^o| + |-^ol- Hence, from ( [T8| ) we notice that finding 
a tight overbounding function Lit{x) as close as possible to p{x), i.e. making \Aq\ as small as possible, 
is crucial for the good performance of a rejection sampling algorithm. 

The most favourable scenario to use the RS algorithm occurs when p{x) is bounded with bounded 
domain. In this case, the proposal PDF, tt{x), can be chosen as a uniform density (the easiest possible 
proposal), and the calculation of the bound L for the ratio p{x)/t:{x) is converted into the problem of 
finding an upper bound for the unnormalized target PDF, p{x), which is in general a much easier task. 
Indeed, in this scenario the performance of the rejection sampler can be easily improved using adaptive 
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schemes |Gilks and Wild 1992 Martino and Miguez 201 la | or strip methods | Hermann et al. 2003 
Chapter 5], [ |Ahrens| [T993| [1995] |Devroye| [T9841 |H6rmann[ [20021 [Marsaglia and Tsang[ [20001 among 
other techniques. Unfortunately, when p{x) is unbounded or its domain is infinite, the proposal 7r{x) 
cannot be a uniform density and, in general, it is not straightforward to design a good proposal PDF (i.e. 
a proposal from which samples can be easily drawn and with a shape as close as possible to the shape of 
the target PDF) inside an infinite domain | Devroye| 1986 Goriir and Teh 201 1[ Hormann et al.[ 2003 [ 



Martino and Mfguez 2011b|. 



Figure [4] illustrates the three possible cases considered in the sequel: bounded PDF with an infinite 
support. Figure [4ja), unbounded PDF with a finite support. Figure [4jb), and bounded PDF with a finite 
support. Figure Qc). In fact, there exists a fourth possible scenario: an unbounded PDF with an infinite 
support. However, since we can consider this case as a combination of the first two cases shown in Figure 
[4j^a) and Figure |4jb), it will only be briefly discussed. The next two sections are devoted to describing 
methods that deal with these problematic situations by transforming p{x) and embedding it inside a finite 



region. First, Section III describes the transformed rejection (TR) sampling approach, and then Section 
IV describes the generalized ratio of uniforms (GRoU) technique. 




(a) 



(b) 



(c) 



Fig. 4. Three possible cases of unnormalized target density, p{x), with three possible overbounding functions Lir{x): (a) 
bounded with infinite domain, (b) unbounded inside a finite domain and (c) bounded with a finite domain. Only in the last case 
it is possible to use a uniform distribution as a proposal PDF, n{x). 



HI. Transformed rejection method 



As already discussed in Section II-D the simplest scenario for the RS algorithm occurs when the 
target density is bounded with bounded support, since a uniform PDF can be used as proposal density, 
as suggested by several authors (see e.g. |Botts et aL| [201 1[ |Devroye[ |1986[ [Hormann] [1993[ [Hormann 



and Derflinger[ [T994 Marsaglia 1984[ Wallace 1976 1). Therefore, an interesting and very active line of 
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research is trying to find a suitable invertible transformation of the target RV that allows us to apply RS 
to a bounded PDF defined inside a finite domain, where we can use the uniform or some other simple 
proposal. Namely, our goal is finding a transformation that converts a PDF of the type displayed in Figure 
Qa) or Figure |4jb) into a PDF of the type depicted in Figure |4jc). 

Conceptually, we can distinguish two cases: a bounded target PDF defined inside an unbounded domain, 
as in Figure Q a), and an unbounded target PDF with bounded support, as in Figure Qb). The third case, 
unbounded target PDF with unbounded support, can be dealt with as a combination of the other two 
cases. Moreover, we can tackle the problem by applying a transformation directly to an RV distributed 
according to the target PDF, X ~ p{x), or to an RV that follows the inverse target PDF, Y ~ p^^{y)- 
Hence, taking into account all the possibilities, in the sequel we have to consider six different situations: 

A. Applying a suitable invertible transformation to an RV X ~ p{x), obtaining Z = f{X) ~ p{z). 

1) When p{x) is bounded with unbounded domain. 

2) When p{x) is unbounded but has a finite support. 

3) When p{x) is unbounded and with an infinite support. 

B. Applying an appropriate invertible transformation to an RV Y ~ p~^{y), obtaining U = h{Y) ~ q{u). 

1) When p{x) is bounded with unbounded domain, implying that p^^{y) is unbounded but with 
bounded support. 

2) When p{x) is unbounded but has a finite support, implying that p^^{y) is bounded but has an 
infinite support. 

3) When both p{x) and p^^{y) are unbounded and with an infinite support. 

Finally, before discussing in detail all these cases, it is important to remark that we can always generate 
samples distributed according to the target PDF from samples of the transformed RVs. On the one 
hand, when an invertible transformation Z = f(X) is applied directly to the target RV, X ~ p{x), the 
transformed RV follows an unnormalized PDF p{z) = p{f~^{z))\f^^{z)\, and, given a sample z' from 
p{z), then x' = f^^{z') is clearly distributed as po{x) oc p{x). On the other hand, if the invertible 
transformation U = h{Y) is applied instead to the inverse target RV, Y ~ p~^{y), then the resulting RV 
U follows an unnormalized PDF q{u) = p^^ {h^^ {u))\h^^ {u)\. Unfortunately, the relationship between 
samples u' from q{u) and samples x' from pq{x) oc p{x) is not trivial, but can still be found and exploited 
to obtain samples from the target PDF, as shown in Section |V] 
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A. Transformation of the target random variable X 

In this section we look for suitable transformations, f{x), applied directly to the target RV, 
X ~ Pq{x) cx p{x), such that the resulting RV, Z = f{X) ~ p{z) is bounded with bounded support. 
In the sequel we will consider, without loss of generality, that f{x) is a class monotonic (either 
increasing or decreasing) function inside the range of interest (i.e. inside the domain of the target RV 
X, Dx)]^ This implies that /(x) is invertible, and its inverse, f~^{z), is also a class monotonic 
function inside the range of interest (the domain of the transformed RV Z, T>z = (0, 1), in this case). 
Finally, regarding the unnormalized target PDF, p{x), we do not make any assumption (e.g. we do not 
require that p(x) is neither monotonic nor continuous) and consider a domain T>x = K for PDFs with 
unbounded support (cases 1 and 3) and Vx = {a, b), with a, 6 G M and a < b, for PDFs with bounded 
support (case 2). 

1) Bounded target PDF p{x) with unbounded support: When the target PDF is bounded with 
unbounded domain and the transformation is applied directly to the target RV X, the sampling technique 



obtained is known in the literature as the transformed rejection method, due to [Botts et aLj |2011[ 
Hormann 1993 Hormann and Derflinger 1994[ Wallace[ [1976 1. However, this approach is also called 



the almost exact inversion method in [ Devroye] |1986[ Chapters 3] and the exact approximation method 
in I Marsaglia} 1984 1, remarking its close relationship with the inversion method |Devroye 1986[ Chapter 
2], as explained later. 

Let p{x) be a bounded density with unbounded support, Vx = IR, and let us consider a class 
monotonic transformation, / : M — (0,1). If X is an RV with unnormalized PDF p{x), then the 
transformed random variable Z = f{X) has an unnormalized density 



p{z) = p{f-\z)) 



dz 



p{f~'{z))\r\z)l for zG(0,l), 



(19) 



where f~^{z) is the inverse function of f{x). Thus, the key idea in |Wallace 1976[ is using an RS 
algorithm to draw samples from p{z) and then generating samples from the target PDF by inverting the 
transformation f{x), i.e. drawing z' ~ p{z) and then taking x' = f~^{z'). By choosing an adequate 
transformation f{x), such that p{z) is also bounded, this strategy allows the proposal PDF, 7r(x), to be 
a uniform density, as in Figure |4jc). 



function f{x) is said to be of class if it is continuously differentiable, i.e. if fix) is continuous, differentiable, and its 
derivative, f(x), is also a continuous function. 
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Obviously, the domain of p{z), T>z = (0,1), is bounded. However, in general the density p{z) can 
be unbounded, i.e. it may have vertical asymptotes, depending on the choice of the transformation f{x). 
Indeed, taking a closer look at ([T9]l we notice that, although the first term ^(/"^(-z)) is bounded (since 
Pq{x) is assumed to be bounded), the second term, is unbounded in general, since 



lim 


df-Hz) 


= lim 


dr\z) 


2-5-0 


dz 


Z-5-1 


dz 



This is due to the fact that f{x) must have horizontal asymptotes, since it is a monotonic continuous 
function that converts the infinite support of p{x), Vx, into a finite domain, Dz = (0, 1). Consequently, 
f^^{z) must have vertical asymptotes at the extreme points of Vz, implying that the limits in ( |20] ) diverge 
to infinity. Figure [5] illustrates this situation, showing an example of a non- monotonic unnormalized target 
PDF with support Vx = and two examples of possible transformations f{x) and f^^{z) (strictly 
increasing and decreasing respectively), where the asymptotes can be clearly appreciated. 




(a) (b) (c) 



Fig. 5. (a) A bounded target PDF p{x) with an unbounded domain Dx ~ R- (b) Two possible examples, monotonically 
increasing (solid line) and monotonically decreasing (dashed line), of the transformation f{x) with horizontal asymptotes at 
a; = and x = 1. (c) The corresponding inverse transformations f^^{z) with vertical asymptotes at 2: = and z — 1. 



Hence, as a conclusion, it is clear from ^19\ and ( [20| ) that the unnormalized density p{z) resulting 
from the transformation f{x) remains bounded when the tails of p{x) decay to zero quickly enough, 
namely, faster than the derivative ^-^ ^J^-* = ^^^^jT^^ diverges when z — z* G {0, 1} (or equivalently, 
when X — )■ /^^(z*) = ±00). More formally, let us note that the limit of interest can be expressed as 

Li = lim,p(z) = limXrHz))\r\z)\ = liiTi^ ^t^^^ = lim (21) 
for z* € {0, 1}, with both ^(/^^(z)) and \ f{x)\^_,i,^^ tending to zero as z — )• z*. Hence, this limit 



will be finite if and only if p(/ ^(z)) is an infinitesimal of the same or higher order than 
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at z = z*. Alternatively, using the last expression of the limit, Li will be finite if and only if p{x) is an 
infinitesimal of the same or higher order than |/(x)| at x = f~^{z*) = ±00. 

We also remark that /^^(z) has vertical asymptotes at both extreme points of Vz because the support 
considered for the target RV X, Vx, is a bi-infinite interval (i.e. it extends towards infinity in both 
directions). If the support of X is a semi-infinite interval in R (i.e. an interval that extends towards 
infinity only in one direction, e.g. Vx = or Vx = IR^), then f^^{z) only has one vertical asymptote 
either at z = or at 2; = 1, depending on the open end of the interval and on whether f{x) is increasing 
or decreasing. However, by focusing on the single asymptote of f~^{z), the discussion performed above 
remains valid. 

Finally, it is also important to realize that better acceptance rates can be obtained by a suitable choice 
of the transformation function /(x). Indeed, when f{x) is similar to the unnormalized CDF, Fx(x), the 
PDF p{z) becomes flatter and closer to a uniform density, so that the acceptance rate using a uniform 
proposal, 7r(z) = U{Vz), is improved. In fact, if f{x) is exactly equal to the unnormalized CDF of X, 
i.e. f{x) = Fx{x), then p{z) is the uniform density inside the interval Vz = [0, 1]. For this reason, this 
technique is also termed almost exact inversion method by some authors (see e.g. | Devroye| 1986 1). 



2) Unbounded target PDF p{x) with bounded support: A similar methodology can also be applied 
when the target PDF, p{x), is unbounded but has a bounded support, Vx = {a, b). In this case, using 
again a class monotonic transformation, / : {a, b) — (0, 1), we can also transform p{x) into a bounded 
density with bounded domain, Vz = (0, 1). For ease of exposition, and without loss of generality, let 
us assume that p{x) has only one vertical asymptote at x = x* , i.e. < p{x) < 00 for x / x* and 
\\m.x^x* fix) = 00, as illustrated in Figure Qb), where x* = a|^Now, let us consider a target RV X 
with PDF p{x) and Z = f{X). We already know that the unnormalized density of Z is given by 

p{z) = p{r\z))\f-\z)l for ze(0,l). (22) 

Unfortunately, although |/~^(z)| is bounded (since f^^{z) is a class function), p{z) is unbounded 
in general, as the first term diverges, i.e. 

lim = lim =00. (23) 

^In many cases, the vertical asymptote of p{x) is located at one of the extreme points of the support. Hence, for a monotonically 
decreasing target PDF with support T>x = (a, b] we will typically have x* = a, as shown in Figure |4|b). 
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Following a similar line of reasoning as in the previous section, we notice that now the limit of interest 
is given by 



p{x) 



La = lim 0(2)= lim p{f ^{z))\f ^iz)\ = lun p{x)\f ^iz)\ = lim 



Thus, since p{x) 



z^f(x*) 

00 when X 



(24) 



X* , a necessary condition for obtaining L2 < 00 is having 



1/ 



-1/ 



\z=f{x) 



0, or equivalently |/(x)| — )• 00, when x ^ x* . However, this condition is not sufficient 



for ensuring that ( |24| ) is bounded. Focusing on the last expression of this limit, we notice that L2 will 
be finite if and only if l/|/(x)| is an infinitesimal of equal or higher order than \/p{x) at x = x* . 

Figure [6] shows an example of an unbounded target PDF, p{x), with a bounded support, Vx = (a, h], as 
well as an adequate transformation f{x) that allows us to achieve a bounded unnormalized transformed 
PDF, p{z). Finally, let us remark that, for a more general unnormalized target PDF, p{x), with several 




(a) 




(b) 



Fig. 6. a) Example of an unbounded target PDF, p{x), with a bounded support Dx = {a, b]. b) Example of a suitable increasing 
transformation f{x) fulfilling that f{x) — > 00 when x ^ x* — a faster than p{x). 



vertical asymptotes located at x* G X*, the same restrictions apply. Indeed, p(z) will be finite if and 
only if l/|/(x)| is an infinitesimal of equal or higher order than l/p{x) at all x = x* ^ X* . 

3) Unbounded target PDF p{x) with unbounded support: Combining the reasoning followed in the 
previous two subsections, it is straightforward to see that the necessary and sufficient conditions for 
obtaining a bounded PDF, p{z), when the target PDF, p{x), is unbounded and has an unbounded support, 
are: 

1) The target PDF, p(3;), must be an infinitesimal of the same or higher order than |/(x)| at x = f~^{z*) 
for all z* G Z* and Z* denoting the set of vertical asymptotes of f~^{z). 

2) must be an infinitesimal of the same or higher order than l/p{x) at x = x* for all x* G X* 
and X* denoting the set of vertical asymptotes of p{x). 
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B. Transformation of the inverse target random variable Y 

In this section we perform the complementary study of the one shown in Section |III-A[ analyzing 
suitable transformations, h{y), applied to the unnormalized inverse target RV, Y ~ p~^{y), such that the 
resulting RV, U = h{Y) q{u) is bounded with bounded support. With respect to the transformations 
we will consider the same restrictions as in the previous section: h{y) belongs to the set of class 
monotonic functions inside the range of interest (i.e. inside the domain of Y , Vy)- Once more, this means 
that h{y) is invertible, and its inverse, /i^^(n), is also a class monotonic function inside the range of 
interest: the domain of the transformed RV U, = (0, 1). Finally, regarding the unnormalized target 
PDF, p{x), now we assume, without loss of generality, that it is monotonic and strictly decreasing with 
a domain Vx = IR^ for target PDFs with unbounded support (cases 1 and 3) and Vx = (0, b] for PDFs 
with bounded support (case 2)|^This ensures that p{x) is invertible and p^^{y) is also a well-defined 
monotonic and strictly decreasing PDF]^ with a domain Vy = for unbounded target PDFs (cases 2 
and 3) and Vy = (0, 1] for bounded target PDFs (case 1)^ 

1) Bounded target PDF p{x) with unbounded support: Let us consider a monotonically decreasing 
and bounded unnormalized target PDF, p{x), with unbounded support, Vx = IR^, such thatp(O) = 1 and 
p{x) — when x — oo. This implies that the unnormalized inverse target PDF, p^^iy), is unbounded, but 
with bounded support, Vy = (0, 1]. Let us consider another RV, U = h{Y) with Y ~ p^^{y), obtained 
applying a monotonic transformation, h{y), bounded inside the domain ofY,h: (0,1] — >■ (0,1). The 
unnormalized density of U is then given by 

g(n) =p-i(/i-i(n))|/i-i(n)| for ne^(0,l). (25) 

Now, since p^^{y) is unbounded when y — y* = 0, in order to obtain a bounded PDF, q{u), a necessary 
condition is 

lim \h'^{u)\= lim |/i(y)r^ = 0. (26) 

u^h(y') J/-5>J/*=0 

''The same conclusions can be obtained using a monotonically decreasing target PDF. However, since monotonically decreasing 
PDFs are more frequently used, we have chosen to work with this class of PDFs. 

^The discussion performed in the sequel can be extended to non-monotonic PDFs. However, in this case we must work with 
the generalized inverse PDF, Pg^(t/), which may be difficult to define in some cases. Thus, for the sake of simplicity we focus 
on monotonic PDFs here, leaving the discussions related to non-monotonic PDFs for Section 



VIII 



Note that using Vy ~ (0, 1] implies assuming that the normalization constant, K, is chosen in such a way that p{0) = 1 
and p(x) — > when x — > oo. 
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However, as it happened in Section III-A2 this is not a sufficient condition. Once more, focusing on the 
limit of interest in this case, 

L3= lim q{u)= lim p-^(h-^{u))\h-^{u)\ = lim P_(^LJ^ = (27) 

u-*h(y) u^Ky) u-^h{y') y^y'=o \h{y)\ 

we realize that a necessary and sufficient condition is that l/|/i(y)| is an infinitesimal of equal or higher 
order than l/p^^{y) at y = y* = 0. 

2) Unbounded target PDF p{x) with bounded support: Here we consider the complementary case 
of the one discussed in the previous section: an unbounded monotonically decreasing unnormalized 
target PDF, p{x), with a vertical asymptote at x = a;* = (i.e. limx^x-'=op{x) = oo), but bounded 
support, Vx = (0,6]. Hence, the unnormalized inverse target PDF, p~^{y), is monotonically decreasing 
and bounded (0 < p^^{y) < b), but has an unbounded support, Vy = M+. Now, considering another 
RV, U = h{Y) with Y ~ P^^{y), obtained applying a continuous and monotonic (either increasing 
or decreasing) transformation, h : IR+ — )■ (0, 1), to the unnormalized inverse target RV, U has an 
unnormalized PDF 

q{u) =p-^{h-^{u))\h-^{u)\ for nG(0,l). (28) 

Again, although the first term, p"^ (^h~^ (u)) , is bounded (as p^^{y) is bounded), this PDF may be 
unbounded, since the second term will be unbounded in general. Indeed, since h{y) transforms an infinite 
domain, Vy = K^, into a finite domain, = (0, 1), it must reach a horizontal asymptote when y — oo. 
This results in a vertical asymptote for /i^^(n) either at u = n* = 1 (when h^^{u) is increasing) or at 
-u = u* = (when h~^{u) is decreasing), implying that 

lim |/i~^(it)| = lim |^(y)| ^ = oo. (29) 



The limit of interest in this case is given by 



L4= lim q{u)= lim p'\h-^{u))\h^^{u)\ = lim = lim , (30) 

u^u. u^u' |/i(y)|^^^_,^_^ y^h-^{u-) \h{y)\ 

Therefore, a necessary and sufficient condition for having L4 < oo is that p~^{y) is an infinitesimal of 
equal or higher order than \h{y) \ at y = h^^{u*). 

3) Unbounded target PDF p{x) with unbounded support: In the more general case {p{x) unbounded 
and with infinite support, Vx = IR^), combining the results obtained in the previous two subsections, it 
can be easily demonstrated that the PDF of the transformed RV U, q{u), will be bounded if and only if: 

1) The unnormalized inverse target PDF, p^^{y), is an infinitesimal of the same or higher order than 
\h{y)\ at y = h^^{u*) for all u* G U* and U* denoting the set of vertical asymptotes of /i^^(tt). 
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2) l/|/i(y)| is an infinitesimal of the same or higher order than \/p ^{y) at y = y* for all y* G 3^* 
and y* denoting the set of vertical asymptotes of p^^{y)- 



C. Summary of the conditions for all the possible cases 

Table |I] summarizes all the possible cases considered in the previous subsections, showing the 
restrictions imposed both on the transformation (/(x) or h{y)) and the target PDF {p{x) or p~^{y)), the 
vertical asymptotes (again both for the transformation and the target PDF) and the conditions required 
for attaining a bounded transformed PDF, either p{z) as given by ( fTQ] ) or p2] ) for Z = f{X) or q{u) as 
given by (125]) or ([28]) for U = h{Y). 



IV. Generalized ratio of uniforms method (GRoU) 



A general version of the standard ratio of uniforms (RoU) method, proposed in |Kinderman and 
Monahan 1977 1, can be established with the following theorem [ Wakefield et al.[[l991[ . 

Theorem 2: Let g{u) : IR+ IR+ be a strictly increasing (in M+\{0} = (0, +oo)) differentiable 
function such that g{0) = and let p{x) > be a PDF known only up to a proportionality constant. 
Assume that {v,u) € M? is a sample drawn from the uniform distribution on the set 



-^0 



(v,n) G : < n < g'^ 



cp\ 



9{u) 



(31) 



where c > is a positive constant and = ■ Then x = is a sample from po{x). 

The proof can be found in the Appendix [a] Choosing g{u) = \v?, we come back to the standard RoU 



method. Other generalizations of the RoU method can be found in the literature [Chung and Lee 1997 



Jones and Lunn[ 1996 1. Moreover, in Appendices [D| and [E| we provide extensions of the GRoU relaxing 
some assumptions. Further developments involving ratio of RV's can be found in ]Barbu[ |1982[ |Curtiss 



T94T] |Dieter| [19891 P^saglia| [T965| [Perez et al.| [2008] [Stefanescu and Vaduval [7987] [Vaduva| [1982 1. 



In literature, the GRoU is also combined with MCMC techniques | Groendykq 2008 1. 



The theorem above provides a way to generate samples from pq{x). Indeed, if we are able to draw 
uniformly a point (u', n') from Ag, then the sample x' = v' /g{u') is distributed according to po{x) oc p{x). 
Therefore, the efficiency of the (standard or generalized) RoU methods depend on the ease with which 
we can generate points uniformly within the region Ag. For this reason, the cases of practical interest 
are those in which the region Ag is bounded. Moreover, observe that if g{u) = u and c = 1 we come 



back to the fundamental theorem of simulation described in Section II-B since Ag becomes exactly ^o- 
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TABLE I 

Summary of the conditions required for attaining a bounded transformed PDF 



conditions 


vertical asymptotes 


conditions tor 


transformation 


PDF 


transformation 


PDF 


bounded PDF 


Z = f{X) 


p(a;) bounded 


\f-\z)\^^ 


None 


p(a;) — > faster than 


/:K^(0,1) 


Ox =K 


when z ^ z* £ {0, 1} 




|/(x)j ^ at X- = / '{z ) 


/ class monotonic 








with z* G {0, 1}. 


z = fix) 


p(x) unbounded 


None 


p{x) — >■ 00 


faster than 


f : {a,b) (0, 1) 


X»x = (a, &> 




when a; — ^ a;* 


l/p(X) —y{)SLtX = X . 


J class ( nionolonic 


a, 6 e R, < /; 








z = f{X) 


p{x) unbounded 


\f-\z)\^oo 


p{x) — oo 


— y faster than 


/ : R ^ (0, 1) 


•Dx =R 


when « — >■ a* 


when a; — > a;* 


|/(a;)|^Oatx = /-i(^*). 


/ class monotonic 








1/1/ I X ) — ? U iaSLCl llldll 

l/p(a^) — )■ at X = X'*. 


U = h(Y) 


p~^{y) monotonically 


None 


p~^{y) oo 


l/|/!-(2y)| — > faster than 


/i: (0,1] ^ (0,1) 


decreasing 




when y — >• y* = 


l/p'\y)^Oaty = y*=0. 


h class monotonic 


Vy = (0, 1] 

unbounded 








U = h{Y) 


p~^{y) monotonically 


|/i~^(m)| -J- oo 


None 


P~^{y) faster than 


h : R+ ^ (0, 1) 


decreasing 


when M — ^ M* 




\h{y)\^Qaty = h-\u*). 


h class monotonic 


T>Y =R+ 
bounded 


with M* = or w* = 1 






U = hiY) 


p~^{y) monotonically 


|/i~"^(m)| — >• 00 


p'^iy) 00 


^/\Hy){ faster than 


h:R+ ^ (0, 1) 


decreasing 


when w — > w* 


when y ^ y* 


l/p~\y)->OsAy = y*. 


h class monotonic 


Vy =r+ 

unbounded 






p~^(j/) — ^ faster than 

/')(//) al // = /( ''^{Jr). 



Note that in the boundary of the region Ag we have u = g ^[cp{x)] and, since v = xg{u), we also 
have V = xg[g~^{cp{x))]. The contour of Ag is described parametrically by the following two equations 

u = g~^[cp{x)], 

V = xg[g~^{cp{x))], 



(32) 



23 



where x plays the role of a parameter. Hence, if the two functions g '^[cp{x)] and xg[g ^{cp{x))] are 
bounded, the region Ag is embedded in the rectangular region 

7^g = |(f,n) e : < n < sup c/"^[cp(x)], 

mi xg[g~'^{cp{x))\ <v< sn^) xg[g^^ {cp{x)y\ \. 

X X ^ 

Figure [7] depicts a generic bounded region Ag, embedded in the rectangular region TZg, defined above. 



(33) 



u 




X 


Kg 


1 ~;::,=-=-^ 






> V 


inixg[g~'^{cp{x))\ 

X 


sup xg[g-^(cp{x))\ 

X 



Fig. 7. Example of a bounded region Ag, embedded in the rectangular region TZg. 

Once, for instance, the rectangle TZg is constructed, it is straightforward to draw uniformly from Ag by 
rejection sampling: simply draw uniformly from TZg and then check whether the candidate point belongs 
to Ag. Note that to use this rejection procedure we do not need to know the analytical expression of 
the boundary of the region Ag, i.e., it is not necessary to know the analytical relationship between the 
variables v and u that describes the contour of the Ag. Indeed, Eq. (31 1 provides a way to check whether 



a point {v, u) e falls inside Ag or not and this is enough to apply a RS scheme. 

Figure jsj^b) provides an example in which the region Ag (obtained with g{u) = and c = 1/2, i.e., 
with the standard RoU method) corresponds to standard Gaussian density (shown in Figure [8ja)). The 
pictures also illustrate different lines corresponding to x constant (dotted line), y constant (dashed line), 
V constant (solid line) in the domain x — y and in the transformed domain v — u. 

In the next section, we obtain the conditions that the function g{u) has to satisfy in order that g^^[cp{x)\ 
and xg[g^^ {cp{x))] be bounded. 



A. Conditions to obtain a bounded Ag 

The region Ag is bounded if the two functions g^^[cp{x)] and xg[g^^{cp{x))] are bounded. Now, we 
study the conditions that the functions g{u) and p{x) have to fulfill in order to obtain u = 5'~^[cp(x)] 
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y 




(a) (b) 

Fig. 8. Examples of the regions Ag obtained applying the standard RoU transformation to standard Gaussian PDF. (a) A 
standard Gaussian density poi^:) oc p{x) — expl—x^ /2}. (b) The region Ag corresponding to a standard Gaussian PDF, 
obtained using g{u) = |it^ and c — I. 



and V = xg[g ^{cp{x))] bounded. 

1) First function u = g^^[cp{x)]: since g is an increasing (g > 0) and continuous function, g^^ is also 
increasing so that the function g^^[cp{x)] is bounded if, and only if, p{x) is bounded, i.e., 

p{x) < M, (34) 

for all G "D, where M is a constant. 

2) Second function v = xg[g~^{cp{x))]: since g > and hence g~^ is also increasing, the function 
xg[g~^{cp{x))] is bounded if: 

a) p{x) is bounded, i.e., p{x) < M, 

b) and the limits 

lim xg[g^^{cp{x))] = Li < +oo, (35) 
lim xg[g~^{cp{x))] = L2 < +00, (36) 



are both finite. The attainment of Eqs. ([35])-([36]) entails the following conditions: 
• Since the first factor in x • g[g^^{cp{x))] is x, we need that g o g^^ o cp vanishes as x — )• ±00, 
i.e., 

lim g[g-\cp{x))] = 0. (37) 
Hence, since lim^—^-i-oo p{x) = (j; G P = M), g ^(0) = (we have assumed ^(O) = in the 



GRoU), and u = g ^[cp{x)] (see Eq. (|32])), the Umit in Eq. (37i becomes 



lim g[g ^{cp{x))] = lim 17 (n) = lim — ^ = 0. (38) 



u-s>0 
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Moreover, since we desire Eqs. (|35])-(|36l), it is also necessary that this factor g[g^^{cp{x))] = 
g{u), for n — )• 0, must decay to zero equal or faster than 1/x — when x — )■ ±00. 
This condition can be rewritten in other forms. For instance, setting y = p{x) (and assuming 
now p{x) invertible, for instance, monotonic decreasing) we can rewrite g[g~^{cp{x))] as 



g[g ^icp{x))] = ^ 



u=g~^{cp{x)) 



dg_ 

du 



.1 . (39) 

u=g-^{cy) 



dg-^ 




dy 


cy 



(recall that c is just a constant) hence when x — )• +00, y = p{x) 0, we need that > 00, 

for y — )• 0, must diverge equal or faster than x = p^^{y) — )■ +00. That is equivalent to assert 
g-lf^y^ — ^ when y — )• equal or faster than 1/x — )• for x — )• +00. Since x = p^^{y) we can 
rewrite it as ^^rr^ vanishes equal or faster than ^pi^ 0, both for y — . 

B. Summary of conditions 

The region Ag generated by GRoU is bounded if: 

1) The function y = p{x) is bounded (i.e., if p{x) is monotonic, x = p^^{y) has finite support). 

2) the limit lim„_5.o ^ = 0, is verified. Since u = g~^[cp{x)] and we set y = p{x), this limit is 



>0 du 

equivalent to limy_>o ^ — = 00, as written in Eq. (39) 



3) The derivative — ^ when n — )• 0, has to vanish to zero equal or faster than 1/x — )• for 



dy 

du ~^ 

X — )• ±00. Setting y = p{x) and consider a monotonic p{x) (so that we can write x = p^^{y)), it 
is equivalent to assert that g-l^y-^ vanishes equal or faster than ^pr^ 0, for y — )■ 0. 
Moreover, we recall that in the GRoU theorem also assumes other conditions over the function g{u): 

4) g{u) must be increasing, 

5) giu) : ]R+ ^ M+, 

6) g{0) = 0. 

We will show that these 3 last conditions can be relaxed. Indeed, they are used to prove the GRoU (see 
Appendix [a]), however they are not conditions needed to obtain a bounded region Ag. 

V. Extended Inverse of density method 



The standard inverse-of-density (loD) method of Section II-C provides the relationship between a RV 



Y distributed as a PDF proportional to p ^(y) and the RV X with a PDF proportional to p(x)|^In 
^Recall that we refer to p{x) and p~^{y) as densities although they are unnormahzed. 
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this section, we study the connection between a transformed random variable U = h{Y), where Y is 
distributed according to p~^{y), and the random variable X with PDF p{x). 

Given a random variable Y with PDF p^^{y) and IJ = h{Y), where /i is a monotonic function, we 
know that the density of IJ is 



dh 



-1 



du 



(40) 



q{u)=p-\h-\u)) 

Denoting as Ah the area below q{u) (see Figure [vjb)), our goal is now to find the relationship between 
the pair ([/, V) uniformly distributed on Ah and the RV X with density p{x). 

It is important to observe that {U, V) = {U, V) since, for the fundamental theorem of simulation (see 
Section II-B i, if {U, V) is uniformly distributed on Ah then U has pdf q{u) so that U = U. Hence, for 



lack of simplicity, in the sequel we use u instead of u, and U instead of U. Obviously, if we are able to 
draw a sample u' from qiu), we can easily generate a sample y' from p^^{y) as 

y' = h-\u'). (41) 

Therefore, using the inverse-of-density method, we can obtain a sample x' from p{x) as 

x' = z'p~\y') = z'p-\h-\u')), (42) 



where z' ~ U{[Q, 1]), u' ~ q{u) and y' = h^^{u') from Eq. (41 1. Equation (42i above connects he RV's 
U and X. However, we are looking for a relationship involving also the random variable V. 

Moreover, we denote as Ah the region delimited by the curve v = q{u) and the axis u. Figure [vj^b) 
illustrates the PDF q{u), the area Ah and a point {u',v') drawn uniformly from Ah- To draw a point 
(u', v') uniformly from Ah, we can first draw a sample u' from q{u) and then v' uniformly the interval 
[0,q{u')], i.e., v' U{[0, q{u')]). Therefore, the sample v' can be also expressed as 



V = zq{u ), 



where z' ~ Z^([0, 1]). Substituting q{u) in Eq. (40l into Eq. (43 1, we obtain 



z'p-\h-\u')) 



dh- 



du 



(43) 



(44) 



q(u') 



Furthermore, recalling Eq. (42 1 we can see that 

v' = z'p^^{h^^{u')) 



dh-^ 



du 



(45) 



u' 



hence 



dh-^ 
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Then, finally we can also write 



I dh-^ 
I du 



v'\h{h-\u'))\, 



(47) 



that is a sample from p{x). We indicate with h = ^ the first derivative of h{x). Eq. (47 1 can be also 



seen as an extension of the fundamental theorem of simulation (Section II-B[ ). 





(J) 










- — ( 































u=h(y) 






(b) 





(a) 

Fig. 9. (a) Given a point (a;', y') uniformly distributed on ^o, y' has PDF p^^{y) wliile x' is distributed as p{x), as affirmed 
by fundamental theorem of simulation and the inverse-of-density-method. (b) Given a transformation of RV U = h{Y) with 
PDF q{u), and a point (u , v') uniformly distributed on the area Ah below q{u), then the sample x' — v' h{h~^ (u)) has density 
p{x). 



Figure |vja) depicts the area delimited by p{x) and a point {y',x') drawn uniformly from ^o- As 



explained in Section II-C x' is distributed as p{x) and y' is distributed as p ^{y). For the standard loD 
method we know that X = Zp-^{Y) where Z ~ ZY([0, 1]) and Y ~ p'^{y). 

Moreover, Equation (47 1 connects a uniform random point (f7, V) € Ah, as illustrated in Figure |vjb), 
and the RV X. Therefore, if we are able to draw points {u',v') uniformly from Ah we can generate 



sample x' from the density p{x) using Eq. (47 1, as formalized by the following proposition. 

Proposition 1: Let y be a RV with a mono tonic PDF p^^{y), and let U = h{Y) be another 
(transformed) RV, where h{y) is a monotonic transformation. Let us denote with q{u) the density of 
U and let Ah be the area below q{u). If we are able to draw a point {u',v') uniformly from the region 
Ah, then 

v' 

^' = = v'\h{h-\u))\, 

I du 



is a sample from the PDF p{x). 
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Below, we provide two interesting special cases. 

• Choosing h{y) = y (hence h = 1), we have U = Y and as a consequence q{u) = q{y) = p^^{y) 



and the region Ah is exactly ^o> so that Eq. (47) becomes 

x' = v', (48) 

i.e., we come back to the fundamental theorem of simulation. Indeed, if we are able to draw a 
point (n' = x',v' = y') uniformly from = -^h, for the fundamental theorem of simulation, it 
yields that x' = u' has PDF p[x) while, clearly y' = v' has distributed as the inverse PDF p-^{y) 
(consideration used in the standard loD method). 
• Moreover, if we take h{y) = \/2y, y > 0, since h~^{u) = ^u^, we have 

x' = ^, (49) 

w 

that corresponds to the standard RoU method. 

VI. Relationship between the GRoU, transformed rejection and IoD methods 
This section is devoted to expose the following proposition. 

Proposition 2: The generalized RoU method can be seen as a combination of the transformed rejection 
method applied to random variable Y distributed according to the inverse density p^^{y), described in 
Section III-A2[ with the extended inverse-of-density method explained in Section |V] 



We first investigate the connection between GRoU and transformed rejection, and then the connection 
between GRoU and the inverse-of-density. 

A. Connection between GRoU and transformed rejection 



Let us recall the region defined by the GRoU in the Eq. (31 1 (for simplicity in the treatment we set 
c = l) 

Ag = \{v,u)eR^ -.Ql^ui^g-^ p[j^ }, (50) 

where, p{x) oc po{x), g{u) is a increasing function and ^(0) = 0. 
I) We consider first for lack of simplicity a monotonic decreasing bounded target density y = p{x) oc 
Pq{x) with an unbounded support V = [0, +oo) (hence the mode is at x = 0). 
Since that g is increasing (then g^^ is also increasing) we can write 
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Moreover, recalling that p{x) = Pdec{x) is decreasing (hence also p ^{y) = Pdec^v) decreasing), 
we have 

and since g{u) > ((7 is increasing), we obtain 

V <p^^^{g{u))g{u). 

Finally, since x S [0,+cxd), then x = p~^{y) > and p~^{g{u))g{u) > (we recall g{u) > 0). 
Therefore, we can write 

0<v<p^^\{g{u))g{u). 

Then these trivial calculations lead us to express the set Ag as 

= {{v,u)eR^ -.OK v<p^^\{g{u))giu)} , (51) 

where p^^^{y) is the inverse of the target density. It is important to remark that the inequalities 
depend on the sign of the first derivative of g (increasing) and p (decreasing). 

2) Similar considerations can be developed for monotonic increasing PDF p{x) = Pinc{x) with an 
unbounded support x G P = (—00, 0] (i.e., x < 0). Indeed, in this case we can rewrite Ag as 

= {(^' n) G m2 . p-^i^ < ^; < 0} , (52) 

Note that since x G (— oo,0], i.e. x < 0, then x = p~^c{y) < and p^^^{g{u)) < then finally 
Pine — 0- The inequalities are different because here p{x) = pinc{x) is increasing. 

3) Similar arguments can be also extended for non-monotonic PDFs. See for instance Figure [8jb) where 
we have p{x) is increasing in Vi = {—00, 0] and p{x) is decreasing in 2?2 = [0, +00), i.e., p{x) is 
non-monotonic with mode located at 0. Moreover, if the mode of the PDF are not located in zero 
o there are several modes, then more but similar considerations are needed. In Section |VIII-C| we 
discuss these more general cases. 

Consider now an increasing differentiable transformation u = h{y) and consider the random 
variable Y with a decreasing PDF p^^{y) and the transformed variable U = h{Y) with density 
q{u) = p-^{h-^{u))h-^{u). The region below q{u) is 

Ah = \{v,u) G : < V <p-^{h-^{u))h-^{u)^ , (53) 



and we can note that Eq. (51 1 is equivalent to Eq. (53 1 when 

y = g{u) = h-Hu). (54) 
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Moreover, clearly, the cases of interest are those in which the region Ag and Ah are bounded, as 
seen in Sections |III-B1| and IV Specifically, in Section III-B 1 we have discussed the properties that 
a transformation h{y) have to fulfill in order to obtain a bounded region Ah, while in Section IV-B we 



have described the conditions to obtain a bounded set Ag. 



It is important to remark that these conditions coincides if we choose g{u) = h ^{u) (with /i(0) = 0, 



see Section |VIII-C| and Appendix |D| about this assumption). Namely, the conditions that the function 
g{u) in Section IV must satisfy in order to guarantee the the region Ag be bounded are exactly the 
same conditions that have to be imposed on the function h^^{u) of Section III-A2 in order to apply the 
transformed rejection method. Therefore, we can state the following result. 

Proposition 3: The region Ag can be obtained as a transformation h = g^^ of a random variable Y 
distributed according to the inverse PDF p~^{y). Specifically, given a RV U = h{Y) = g^^iY) with 
PDF indicated as q{u), the region Ag coincides with the area Ah below the curve q{u). 
This proposition means that the set Ag defined by Eq. ( 3 1 1 or (BTl) is obtained by applying the transformed 



rejection idea for unbounded PDF's to the inverse density p~^{y) (see Section III-A2i. Figure 10 ^b) 
displays the region Ah (that coincides with Ag if g = h"^) defined in Eq. (53 1. Figure 10 c) depicts the 
same region Ah rotated 90°. Proposition [3] can also be deduced as shown in Appendix [B] 




v' 




X - 1 ,1 




1 <'" 1 






/ \q{u) 




Ak \ 









(b) 




(c) 



Fig. 10. (a) Example of region defined by the inverse density p~^{y). (b) The density q{u) = | | p^^ {h~^ (u)) 
obtained transforming the RV Y, i.e., U = h{Y). Generating uniformly the point (u', v') in the area Ah we can obtain samples 
x' from p{x) using Eq. \55\. (c) The region Ah rotated 90° in order to show it how appears when we apply the GRoU technique. 



Moreover, Proposition |3] yields the following corollary about the two marginal densities of the random 
variables U and V with uniform joint pdf on the region Ag provided by the GRoU. 

Corollary 1: Consider a random vector {V, U) uniformly distributed on the region Ag provided by 
the GRoU. We already know that the RV X = -Xj^ has pdf p{x) as proven by the GRoU. Moreover, we 
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can assert that U is distributed as q{u) = p ^{diu))^ (since U = g ^{Y)) and V is distributed as the 



generalized inverse density q ^{v) of q{u) (see Section VIII for the definition of the generaUzed inverse 
pdf). 

B. Connection among GRoU, extended loD and fundamental theorem 

Moreover, In Section |v] we have analyzed the relationship between X with PDF p{x) and the RV 
U = h(Y), where Y is distributed as p^^{y). Hence, given two samples v' and u' uniformly distributed 
on the set Ah, the area below the PDF q{u), we can assert that the sample 



, _ v' _ v' 

is distributed as pq{x) oc p{x), as we prove in Section |v] for the extended loD and extended fundamental 
theorem of simulation. Note that, if we set h{y) = g~^{y), we obtain x' = v'/g{u') that is exactly 
equivalent to the GRoU technique in Section |IV] Therefore, we can also assert the following two 
propositions. 

Proposition 4: The GRoU extends the underlying idea of the classical inverse-of-density approach, 
described in Section |II-C[ Indeed, the classical loD method uses a random variable Y distributed as 
the inverse density p~^{y) to draw samples from p{x), whereas the GRoU uses a transformation of the 
random variable Y, U = g^^{Y), to generate samples from p{x). 

Proposition 5: The GRoU can be also seen as an extension of the fundamental theorem of simulation, 



described in Section II-B Indeed, the fundamental theorem links the coordinates of a random point 



{X,Y) € ^0 with the PDFs p{x) and p ^(y), i.e., X ~ p{x), Y ^ p ^{y), whereas GRoU links 

V 
9(U) 



the coordinates a random point {U,V) G Ag with the same PDFs p{x) and p ^(y), i.e., -^Xjr ~ p{x), 



g{U)^p-\y). 

Clearly, Propositions [3} |4] and |5] entail the Proposition [2] Note that the GRoU can also be extended 
for a decreasing g{u), as we will show in Appendix |e] (see also Appendix |C]l. In this case, we have 

V 

9(u) 



X = then in general for the GRoU we can write 



V 

X 



\9{u)\ 



exactly as in Eq. (55 1. 
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C. Function g{u) to obtain a rectangular region Ag and first formulation of the loD 

Clearly, the easiest case to perform exact sampling with GRoU is that Ag be a rectangular regiorj^ 
The considerations in Section |Vl] are very useful to clarify which g{u) produces a rectangular region 
Ag. More specifically, Proposition |2] allows us to infer which is the optimal (theoretical) choice of the 
function g{u). 

Indeed, since the GRoU is a transformation of a RV Y with pdf p~^iy) (considering, for instance. 



a decreasing p{x)), specifically U = g ^{Y) with c = 1, the well-known inversion method |Devroye 



19861 asserts that if the function g ^{y) is the cumulative distribution function (CDF^ of Y then the 
transformation produces a uniform RV U. Hence, the set Ag is a rectangular region if we use 

g-\y) = FY{y)^g{u) = Fy\u), (56) 
where Fy(y) is the CDF of RV Y, i.e.. 



Fviy) 



p ^{y)dy. 



(57) 



Since p ^{y) is unnormalized, note that Fyiy) 1/K with y — )• +oo (instead of Fyiy) — s- 1) where 

1 



I p ^{y)dy = / pix)dx. 

JVy JVx 

Therefore, if g'^iy) = Fyiy) then U = Fy{Y) is a uniform RV in [0, 1/K], and Ag is a rectangle 
0<u<^andO<i)<l as we show below in Eq. i60[. Indeed, when g ^{y) = Fy{y), the region 
Ag is defined (c = 1) as 



An 



(u, n) G : < u < Fy 



P 



Fy\u) 



and since 



Y 



GR-^ :0<u<Fy 



p{vp-HFyHu))')]} , 



(58) 



(59) 



and X = vp~^{Fy^{u)) is distributed as po{x) oc p{x). Since < Fy{y) < 1/K, then the suitable values 
of the variable u are contained in [0, The variable v is contained in [0,1] independently of the 
values of u, because inverting the inequalities in Eq. ( [59] ) we obtain 

P~HFy\u)) _ ^ 



0<v< 



p-HFyHu)) 



Clearly it is just one possibility, there are other situations where we can perform exact sampling (for instance, if Ag is a 
circle or a triangle). 

'The CDF Fy (y) of RV Y can be easily expressed as function of Fx (x) (the CDF of X) for monotonic decreasing target 
pdfs po{x) oc p{x) as we show in Appendix jpj 
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SO that Ag is completely described by the inequalities 



An 



{v,u) : <v <1, < n < 



K 



(60) 

obtained with 



Moreover, observe that the expression x = vp~^{Fy^{u)) is exactly the same of Eq. (11 
the first formulation of loD method. Indeed, in this case Ag is a rectangle {v,u) G [0,1] x [0, 1/i^] 
then V ~ U{[0, 1]) and the sample y = g{u) = Fy ^(n) is distributed as p^^{y), so that x = vp^^{y) is 
equivalent to Eq. ( [TT] ). 

D. The second formulation of the loD (Khintchine's theorem) as special case of the GRoU 

Here we show that the second formulation of the inverse of density method described at the end of 
Section |II-C| is contained by the generalized ratio of uniforms technique. Assuming an increasing target 
PDF p{x), if we set g{u) = p{u), namely we use as function g exactly our target PDF p, and c = 1 then 



Ag=p = {{v,u) gR'' -.0 <u<p 



-1 



p 



p{u) 



A, 



{v,u) eR"^ ■.0<u< 



p{u) 



(61) 
(62) 



or using the alternative definition of the area Ag in Eq. (51 1, we have 

Ag=p = {{v,u) -.Q <v <up[u)] . (63) 
Hence, the region Ag=p represents the area below the vertical density (see Eq. 13 1 |Jones 2002 



Khintchine 1938 Troutt et al. 2004 1 corresponding to p ^{y) used in the second formulation of the 
inverse of density method (see Section |II-C[ ). The ratio of uniforms approach assert that 

(64) 



p{u') 



is a sample from p{x) if u' is distributed as the vertical density up{u) and v' ~ U{[Q,u'p{u')]). Now, 
we want to certify if this statement is also true using the second formulation of the inverse of density 
technique. 

Consider a sample z' ~ ^([0, 1]) then z'u'p{u') ~ U{[Q,u'p{u')]), hence we can write 



v' = zup{u). 



(65) 



Replacing the relationship above in Eq. ( 64 1 we obtain 



x' = z'w', 



(66) 
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where z' ~ ^([0,1]) and w' is drawn from up{u). Note that Equations (14i and (66 1 coincide since, 



up{u) is exactly the vertical density of p ^{y) in Eq. (13 1. Therefore, we can assert that the second 



formulation of the standard inverse of density in Section II-C can be found choosing g{u) = p{u) in the 
GRoU (where p{u) is our target PDF). 

Figure [TT] summarizes the relationships among densities, random variables and sampling methods (the 
two versions of the loD, the VDR and the GRoU) for a decreasing target PDF p{x). 



Y ~ p-\y)\ '[ElZ| |t/ ~t/([0,/^])| 
z~t/ao,i])» 




if/iO') = /'"'{>' 
J 



U-q{u) = p-\h-'{u))Y^ 



Fig. 11. Relationships among densities, random variables and sampling methods for a decreasing target PDF y = p{x). Dashed 
line shows the connection produced by the GRoU algorithm whereas dotted line describes the second formulation of the inverse 



of density (loD) method. The random variable W ^ w-£^ is, the vertical density corresponding to the pdf (y). Finally, 
dashed-dotted line depicts the relationship between the first formulation of the loD and GRoU when h{y) — g^^{y) = 
(see Section [VTC}. 



E. Effect of the constant c 

So far, for simplicity we have set c = 1. However, all the previous considerations and remarks remain 
valid. Indeed, assuming a decreasing p{x) = Pdec{x) and c > 0, for instance Eq. ( [5T] ) becomes 



= \{v,u)(^M? -.QKvK p-^l ( ^ ) g{u) \ . (67) 
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Since, as we show in Section [TTj all this techniques work with unnormalized PDF, we can also multiply 
both inequalities for a positive constant 1/c obtaining 



^ c I c 

-1 ( <M\ sH 



(68) 



then Ag^p^^^ represents the area below q{u) = p^^^ \c:) c; ^^^^ (unnormalized) PDF of the RV 
U = g~^{cY) where Y is distributed according to (unnormalized) PDF p^^f.{y)- See also Appendix |b[ 



VII. Further Considerations 

In this section, we provide further observations about the connection among the GRoU and the 
transformed rejection sampling, and about some assumptions over the function g{u). 

A. Minimal rectangular region 

Here, we show that the minimal rectangle TZg such that Ag C TZg is equivalent to a rectangular region 
TZh embedding a set Ah obtained with a trasformation of a random variable Y with density p^^{y). 

We have seen that the minimal rectangular region embedding the region Ag {Ag C TZg) of the GRoU 
is defined as 

7^g = |(7;,n) G : < n < sup c/"^[cp(x)], 

(69) 

m.{ xg[g'''^{cp{x))\ <v< s\i]ixg[g~^{cp{x))\ \. 

X X ^ 

For lack of simplicity, in the following we consider a bounded decreasing PDF p{x) defined for all 
X G M+ with mode localized at = 0. Recalling also that g is a positive increasing function, the first 
important observation is that xg[g^^{cp{x))] is also positive so that 

inf xg[g~^{cp{x))] = 0. 

a;eR+ 

Moreover, since g is increasing then also g^^ is increasing, a second observation is that 

sup g''^{cp{x)) = g'^{cp{0)), 

where x = is location of the mode of p{x) (namely p(0) = sup2,g][5+ p{x)). Therefore for a bounded 
p{x) defined in and a mode at 0, we can rewrite TZg as 

7^g = G : < n < sup g~^[cp{0)], 

(70) 

< V < sup xg[g~^ {cp{x))] >. 
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Now, let us consider a increasing transformation h{y) and a random variable Y distributed according 
to the inverse PDF p^^{y). Note that, since we assume a p{x) is decreasing, bounded with mode at 
X = 0, p^^{y) has bounded domain < y < p(0) but it is unbounded with a vertical asymptote at 0. 
Then, we consider the RV U = h{Y) with PDF 

q{u) = p-\h-\u)) , \ with h{0) <u< hip{0)) 
du 

and indicate with Ah the area below q{u). We also assume that h{y) is chosen adequately such that Ah 
is bounded. In this case, a minimal rectangle TZh embedding Ah exists and clearly it is 

TZh = I (v, n) G : h{0) <u< h{p{0)), 

^ . (71) 

(Ih-^fiA I 
0<v< supp'^{h-^{u))- 



du I 

Now, we desire to express q{u) first as a function of y, obtaining q{y), and later as a function of x, 
obtaining q{x). Recall that Y a RV with PDF p^^{y) and U = h{Y). Then, we have u = h{y) and 
y = h^^{u) and we can write 

q{y) =P~Hy)^^^-^^^, with 0<y<p{0). 

Moreover, since x = p^^{y) and y = p{x) we can also write 

, , dh^^(h(n(x))) . , 1 , , 

q{x) = x ^VFW/y^ ^jjj^ < X <p"^(0) ^ +00, 

du 

i.e., 

q{x) = X ■ h~^[h{p{x))], with < x < +oo. 

Then, we can rewrite the minimal rectangle TZh as 

TZh = {iv,u) G : h{0) <u< h{p{0)), 

<v < supx • h^^[h{p{x))], >. 

X J 

and if we choose h = g^^ with /i(0) = and c = 1 then 

TZh = [{v,u) G . Q<u< g-\pm, 

< v < supx • g[g^^ {p{x))]^ \. 



(72) 



(73) 



Note that TZh in Eq. (73 1 is exactly the same rectangle TZg in Eq. (|70|) when a = and c = 1. 
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B. About the condition g{Qi) = 

One assumption of the GRoU is that g{{)) = 0: is it strictly necessary? can this condition be relaxed? 
We can disclose that the condition g{{)) = is needed with the version of the GRoU that we have tackled 



so far, in Section IV However, note that it is possible to propose different versions GRoU as we show 
in the Appendices |Cj |D] and |E] 

To relax the assumption g{{)) = we can study two different possibilities: g{Q) = c 7^ and g{h) = 0. 



The condition (7(0) = c 7^ 0, depicted in Figure 12 a), is impossible (at least, in the classical formulation 



of the GRoU of Section [rV[ ). Indeed, we have that 

• in the standard GRoU the function g{u) has to be increasing, 

• and we know that the transformation h{y) = g~^{y) is applied to a RV Y with PDF p^^{y) that is 
defined in (0,p(0)]. 



In this situation, the inverse function h{y) = g ^(y) is shown in Figure 12 b). We know that this 
transformation U = h{Y) has to be apply to a RV Y with PDF p^^{y) of type in Figure 12 c). Then, 



the RV Y takes values in (0,p(0)]. Therefore, if c 7^ in the interval (0,c] h{y) is not defined but Y 
can take values there, then the transformation U = h{Y) = g^^{Y) is not possible. 




(a) 



(b) 



(c) 



Fig. 12. Explication why g{0) = c 7^ is not possible in the GRoU. (a) Example of (increasing) function g{u) with 
g{0) = c 7^ 0. (b) The inverse function g^^{y) ~ h{y) corresponding to the g{u) in Figure (a), (c) Example of inverse PDF 



However, the second case g{b) = with 6 7^ is possible with a slight extension of the GRoU that we 
show in the Appendix [d| Indeed, in this case we have g{u) : [b, +00) — )■ and then /i(0) = ^^^(O) = b. 
Figure 13 ^a) illustrates an example of function g{u) with g{b) = and b ^ 0. Figure [Tsjb) shows the 



corresponding inverse function g (y) = h{y). Finally, Figure 13 c) depicts the area Ag this case when 
g{b) = then h{0) = b. 
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(a) (b) (c) 

Fig. 13. Explication why g{b) — with & 7^ is possible in an extended version of the GRoU (see Appendix|Dj. (a) Example 
of (increasing) function g{u) with g{b) — with 6 7^ 0. (b) The inverse function g^^iy) = h(ii) corresponding to the g{u) in 
Figure (a), (c) Example of region Ag when g{b) — with 6 7^ 0. 



Another way to understand this issue is the following: in the definition of Ag we need to combine 
g^^ and cp{x), i.e., 

o cp{x) = g''^[cp{x)]. 

Since p{x) G [0,M] where M = snp^fzs^p{x) (M = p{0), we are assuming the mode is localized at 0) 
then g^^ has to be defined in [0, M]. 

C. Constant c and image of g{u) 

The sign of the constant c is related to the domain of g^^{y), namely, the image of g{u). Indeed, if 
c > 0, it is straightforward to see the g^^{y) must be defined in IR+ since we have the composition of 
functions g^^ o cp{x) where p{x) > 0. Indeed, since we have assumed c > 0, so far we have considered 
always functions g{u) : M — )■ M^, i.e., the domain of g~^ is IR+. 



Moreover, as we have seen in Section VI-E for general values of c, the GRoU is equivalent to the 
transformation of RVs U = g^^{cY) where Y has PDF p'^{y), then another time we can deduce that the 
RV cY must take values into the domain of g^^{y). Therefore, if we consider a function g{u) : M — )• M~ 
then we must use a negative c, i.e., c < 0. 



VIII. General PDFs 

In this section, we investigate the connection between GRoU and Inverse-of-density ( Khintchine's 
theorem) for generic densities. 
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A. Inverse-of-density (and Khintchine's theorem) for generic PDFs 

Before to analyze the GRoU applied for generic non-monotonic PDFs, first of all we discuss and recall 



how it is possible apply the inverse-of-density approach in Section II-C for generic PDFs. Let us define 
the set of points 

■^o\y = {{x, z) G ^0, z = y}, (74) 



i.e., all the points in Ao such that 2 = y for all y G 
PDF as 

where |^o|yl is the Lebesgue measure of AQ\y. 



. Then we can define the generalized inverse 

(75) 



Then the inverse-of-density approach (and all the extended versions of Khintchine's theorem [Bryson 
and Johnson j [19821 |Devroye| [T9841 |de Silva| [T9781 [Chaubey et aL] [20T0l [Olshen and Savage| [1970] 



Shepp 19621) can be summarized in this way: we can draw samples from p{x) if we are able 

• to generate a sample y' from p^{y), 

• and then generate uniformly a point {x',y') on ^o|y'- Then x' is distributed according to p{x). 
Note that this approach is strictly related to the slice sampling algorithm. Clearly, this general approach 
can be expressed in different ways in different specific cases (as symmetric unimodal PDF with mode 
at | Shepp[ [1962[ ), yielding different versions of Khintchine's theorem |Chaubey et al. 2010[ Shepp] 
T9621 . 



It is interesting to observe that (a) is monotone non-increasing [Damien and Walker 2001 Jones[[2002| , 
(b) it has an vertical asymptote at (if the domain of p{x) is unbounded) and minimum at sup^p{x). 



Figure 14 shows an example of bimodal PDF and the corresponding generalized inverse PDF pQ^{y). 
Observe that, for instance, in case the set Ao\y can be formed by two disjoint segments (as in Figure 
[I4ja), Si and ^2) or just one depending on the value of y. Clearly, the length of the sets Si and 5*2 
depend on the 4 monotonic pieces Pi{x), i = 1, 4, that form p{x). 



B. GRoU for unimodal PDF with mode at 



In Section VI we have already seen the definition of Ag when p{x) is increasing or decreasing with 
mode at 0. If the target PDF p{x) is unimodal with mode at 0, we can divide the domain V = V1UV2 
where Vi = [0, +cxd) with p{x) = Pdec{x) is decreasing and 2?2 = (— oo,0] with p{x) = Pinc{x) is 



increasing. Hence, in this case, the complete set Ag can be also written as (combining Eq. 51 and Eq 
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>' 'A|,,'=.^^iU5, fx 




(a) (b) 

Fig. 14. (a) A bimodal PDF p(x). Monotonic parts of p{x) are indicated by Pi{x), i — 1,...,4. (b) The corresponding 
generalized inverse PDF pQ^{y). 



52) 

A, = {{v, n) e pr^, {g{u)) g{u) < v < p^^^ {g{u)) g{u)} . (76) 

Note that the inequaUties depend on the sign of the first derivative of p{x) (i.e., where p{x) is increasing 
or decreasing). Then, we could interpret it as if the GRoU applies a transformation U = g^^{Y) over 
two random variables, Yi with PDF p^^^{y) and Y2 with PDF —plnciv)- 

Figure 15 a) shows an example of unimodal PDF with mode localized at zero. Figure [TSjb) illustrates 



the same region ^0 rotated 90° (i.e., switching the axes x and y). It is possible to figure out that 
Pc^iy) — Pdec(y) ~ Pinciy) where PQ^{y) is the generalized inverse density associated to p{x). Finally, 



Figure 15 c) depicts the corresponding region Ag using g{u) = u /2. We can also observe that for a 
given value u and defining the subset Ag\u C Ag 

■^g\u- {{v,z) ^ Ag,Z = u}, (77) 

then we can write the expression 

du 

where in the first side we have the PDF of a transformed RV U = g^^{Y) where Y is distributed as 
PQ^{y) and \Ag^u\ is the Lebesgue measure of the subset Ag^y_- 

C. Unimodal PDF with mode at a ^ 
In this section, we consider the application of GRoU method to a unimodal PDF p{x) (x G IR+ without 



PGi9in))^ = \Agiu\, (78) 



loss of generality) with mode at a / 0. We can see an example in Figure 16 ^a). 
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I 



(a) 



\pl(y) ! 



du 



plliy) 



du 



1 . du 



(b) 



(c) 



Fig. 15. (a) A unimodal PDF pix) with mode localized at zero, (b) The region rotated 90°. In this case, 
the generalized inverse PDF PQ^{y) can be written as PQ^{y) = Pdec(y) ^ PTnciv)- If '^^e target PDF p{x) is 
also symmetric we have p~il^{y) = ^PTnciv)' hence PQ^{y) = 2p'^^^{y). (c) The corresponding region Ag = 
{{v, u)eR^ : pT-„\ (g{u)) g{u) <v< p^^^ {g{u)) g{u)} obtained by the GRoU using g{u) = jl. 



In Figure 16 ^b) is depicted the region below p{x) with the axis x — y switched (as rotated 90°). 



In this case the region can be described as 

A = {(x, y) e M2 : (y) <x< vtM)- (79) 

Moreover, observing Figure [T6jb) we can individuate and define 5 random variables: a RV Y\ with PDF 
Pdec(y) ~ ^ (associated to the region Ai), RV Y2 with PDF a — p^nciv) (associated to the region A2), I3 
with PDF p^nc(y) (associated to the region A2), I4 with PDF p^^^iy) (associated to the regions Ai, A2 
and A2) and finally Ys with PDF the generalized inverse density p^{y) = Pdeciv) ~ Pinciv)- Note that 
Aq is only composed by Ai and A2, i.e., Aq = AiU Ap"| 

Now, we consider the transformation of random variables Ui = g~^{Y^) and U2 = g~^{Yi) (with 
an increasing function) and plot together the two PDFs qi{u) oc Plncidi'^))^ ^"^^ Q2{u) oc P^ec(5(^))af 
obtaining the regions Bi, B2 and as represented in Figure [T6l^c). The region attained with the GRoU 
method is exactly Ag = BiVJ B2 (in Figure 16 c) we use g{u) = n^/2) and we can write it as 



A = [m e : P^\{9{u))f^ < V < P,-c(5(n))^} • (80) 

'"if we desire to draw uniformly on Aq defined as in Figure |l6[b), we should to be able to simulate a RV Y5 with PDF 
Pa^iy)- Indeed, to do it we could simulate a r.v. Y4,, i.e., generate a sample y' according to a PDF proportional to p'^^^iy), then 
draw u ^ W([0, 1]), finally calculate x' — u'p'^^^{y') and accept x' if x' > pJnciv') (hence, {x',y') is uniformly distributed 
on x' is distributed according p{x) and y' as pQ^{y)). 
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Note that we can interpret that the boundary of Ag can be obtained through a "transformation" of the 



contour of the region (see Eqs. (79 1 and (80l). Finally, recalling the subset Ag^u = {{"v, z) & Ag, z 
u} then note that in this case we also have PG^{g{u))^ = 



Pa(y')=pl(y')-p-l(y') 




X 


Paiy') 

A 

1 

1 


A 

a 


1 

At,' -. 




^ y. 



(a) 



(b) 



P<7 («(«)) V 
du 
































u 

• >■ 



du 



(c) 



(d) 



Fig. 16. (a) An example of unimodal density p(x). (b) The region represented switching the axes x — y. (c) The region 
Ag obtained with the GRoU technique with g{u) = (d) The same region Ag represented switching the axes m — w in the 

previous picture (it is typical representation of the GRoU regions). 



D. Generic PDF 

Let us assume that we can divide the domain Vx of the PDF p{x) with a partition formed by A'' 
disjoint sets, i.e., Vx = "Di U 1^2 U ...Vx, where p{x) is monotonic increasing or decreasing, i.e., 

p{x) = pj{x) and xeVj, (81) 

where Pj{x) is an increasing or decreasing function. 

Let us assume, moreover, that p{x) is a continuous function with Vx = Since J^^ p{x)dx < +oo, 
then N is even and ^21-1(2;), with i = l,...,N/2, are increasing functions whereas P2i{x), with 
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i = 1, N/2, are decreasing functions. Then, the region Ag generated by the GRoU can be expressed 
as 

Ag=Ag,lUAg,2^--.^Ag^N/2, (82) 

where 

Ag,i = {iv,u)eM'' : {g{u)) g{u) <v< p"/ {g{u)) g{u) } , (83) 
for 1 = 1, ...,N/2. Figure MTra) shows the bimodal PDF po{x) oc p{x) = exp{ — (x^ — 4)^/4} and the 



corresponding region Ag obtained by the GRoU with g{u) = is illustrated in Figure 17 b). We recall 



that, as illustrated in Figure 17 a), we can define 

A|2/ = {{x,z) eAo,z = y], 

and then we can write 

Since in this case it is composed by two segments, AQ\y = U ^2, we have p^{y) = \Si\ + \S2\ - Then, 
recalling the definition of the subset Ag\u = {(^) ^) ^ -^g^ ^ = hence note that we have again that 

dg 



(84) 



(85) 



du 



as depicted in Figure [TTj^b). 



I \ 
Pi Ml \pii^'} 



P4W 



— I.A 1 1 




u' 








'■-^fi.'t-. 
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^ 



(a) 



(b) 



Fig. 17. (a) A bimodal density po{x) oc p{x) = exp{— (a::^ — 4)^/4} formed by 4 monotonic pieces Pi{x), i — 1, 4. (b) 
The the corresponding region Ag = Ag.i U Ag.2 obtained by the GRoU using g{u) = /2. 
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E. Discussion about GRoU 



The RoU and GRoU techniques were introduced |Kinderman and Monahan 1977 Wakefield et al. 



1991 1 as a bivariate transformation of the bidimensional region below the target PDF p{x). To be 
specific, the RoU techniques were presented as a transformation of a bidimensional uniform random 
variable defined over ^o- This bivariate transformation follows the equations x = and y = u (see 
Appendix [a]). These relationships describe all the points within the transformed region Ag. 

In this work (and, specially, in this section) we have also seen that the GRoU can be interpreted as 
transformations of random variables Yi with PDFs the monotonic pieces p^^{y), i = '^,---,N (where the 
monotonic functions Pi{x), i = 1,..,A^ compose the target density p{x)). These transformed densities 
describe disjoint parts of the boundary of region Ag obtained with the GRoU. 

Furthermore, given the random vector {V, U) uniformly distributed on Ag, we have seen that the second 
random coordinate U is distributed according to q{u) oc {g{u))^. Namely, we can write the RV U 
as a transformation of a RV Y, i.e., exactly as f7 = g^^iY), where Y is distributed according to the 
generalized inverse PDF p^{y). 

IX. GRoU FOR UNBOUNDED PDFS 

Another assumption used on the Theorem |2] of the GRoU is that pq{x) oc p{x) must be bounded. In 
this section, we discuss as to design a GRoU technique for unbounded PDFs (with bounded support, for 
simplicity) using the observations in Section |Vl] We will refer to this technique as unbounded GRoU 
(U-GRoU). Then consider, for instance, a decreasing target PDF pq{x) oc p{x), where 

p{x):Vx = (0,6] 

with an vertical asymptote at x* = 0. In this case, to apply a kind of GRoU approach to draw samples 
from pq{x) oc p{x), we have two possibilities: 



1) the first option is to apply the standard GRoU for bound PDFs of Section IV to the inverse PDF 
p^^{y) (that is clearly bounded with unbounded domain, in this case), in order to produce a sample 
y' from p^^{y). Then, samples distributed according to poix) can be obtained using the loD method 



in Section II-C i.e., x' = z'y' where z' ~ 1]). However, we need to be able to evaluate p^^{y), 
namely to invert p{x), and it could be difficult or impossible, in general. 
2) A more general approach is to design a GRoU technique to tackle directly this kind of unbounded 
target PDFs. To do that, we can use the observations and discussions about the GRoU provided in 



the previous Section VI 
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In Section VI we have emphasized that the GRoU is equivalent to a transformation g {y) of a RV 
Y ~ p^^iy) (with c = 1 and p{x) monotonic). If the transformation U = g^^{Y) is adequately chosen 
the region Ag defined by the GRoU is bounded. 

In this situation, p^^{y) : M+ — )• (0,6] is bounded with unbounded support. In Section III-B2 
have described the conditions that an increasing transformation (p{y) : 



we 



— ^ [di,d2) (where di < d2 
are generic constant) has to fulfill in order to obtain bounded PDFs with bounded support. The random 
variable U = ip{Y) has PDF 



l{u)=p {(f (u)) 



du 



with di < u < d2- 



(86) 



The PDF q{u) is bounded if p~^{y) is is an infinitesimal of the same or higher order than ip{y) at 
y — )• +00, as we have shown in Section III-B2 Therefore, with this suitable function (p{y) and the 
observations in Section VI we can define the corresponding suitable region as 



0<u<ip 



P 



(87) 



so that the sample 



is distributed as po{x) oc p{x) if {v,u) are uniformly distributed on A^. In the sequel, we provide two 
examples of suitable transformations ip{y). 

Example 3: Consider the unbounded target pdf 



Po 



(x) oc p{x) = 21og(x) with x £ (0, 1]. 



(88) 



A first U-GRoU scheme in this case can be found using ip{y) = arctan(y), i.e., if {v,u) are uniformly 



distributed on A^, 



A.II 



< n < arctan 



(89) 



^tan(ii)2 + 1 

then X = tan(M)^+i distributed as po{x). A second possibility is given using, for instance, ip{y) = 
i.e., if {v, u) are uniformly distributed on A^ 



A. 



■v 



{v,u) G 



< li < 



then X 



P \ u/{u~iy^-i/{u-i) 
u/{u-i/-i/(u-i) distributed according to po{x). 



+ 1 



(90) 
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X. Conclusions 



In this work, we have shown that the generaUzed ratio of uniforms (GRoU) algorithm | Kinderman and 



Monahan 1977 Wakefield et al.[ 1991 1 can be seen as combination of other classical sampling strategies: 



an extension of inverse-of-density (loD) method | |Bryson and Johnson 1982[ Chaubey et al. 2010 



Devroyel |1986[ |Jones[ |2002[ |Khintchine| |1938[ |Shepp| |1962[ (and, at the same time, of the fundamental 



theorem of simulation approach |Devroye 1986[ Robert and Casella 2004 1) and the transformed rejection 
sampling (TRS) pevroye| [T9861 [Hormann and Derflinger[ [T9941 |Marsaglia[ [T9841 |Wallace| [T976 |. 
Specifically, for monotonic densities with mode at zero, the GRoU is completely equivalent to a 



combination of the TRS approach applied to the inverse PDF p^^{y) (Section III-Bli, an extension of 
the loD strategy (Section II-C| ) and an extension of the fundamental theorem idea (Section |II-B i. The 
classical loD technique and also the fundamental theorem of simulation provide the relationship between 
random variates from p^^{y) and samples drawn from p{x) whereas the GRoU links the realizations of a 
RV U = g^^{Y), where Y has PDF p^^{y), to the samples distributed as the target PDF p{x). Moreover, 



we have exhibited that the conditions about the function g{u) required in the GRoU [Wakefield et al. 



1991| to obtain a bounded region Ag are exactly the same assumptions needed to the transformation 
U = h(Y) = g^^{Y) (where Y has density p^^{y)) in the TRS approach, in order to achieve bounded 
transformed PDF q{u) = p{h^^{u))^^K-^ with bounded domain. We have also seen that the TRS approach 



can be considered lightly more general than the GRoU approach in the sense that some conditions required 
by the GRoU can be relaxed as we show, for instance, in Sections VII-B||VII-C| and the Appendices D]|E 



For generic non-monotonic densities, formed by N monotonic pieces pi{x), i = I, .., N, the boundary 
of the region Ag of the GRoU can be expressed as transformations of random variables Yi with 
PDFs p^^{y), the inverse functions of the monotonic pieces pi{x), i = 1,..,A^. Moreover, the RV 
U of the random vector {V, U) uniformly distributed on Ag is distributed according to the PDF 
q{u) oc Pci'idi'^))^ where p^{y) is the generalized inverse PDF. Namely, we can write the RV U 
as a transformation of a RV Y, i.e., exactly d&U = g^^{Y), where Y is distributed according to p^^{y). 

Therefore, in this work we have illustrated the close relationships among GRoU, loD and TRS 



approaches in Section |Vl] Using the considerations in Section |Vl} we have also relaxed different 
assumptions of the GRoU (see Sections VII-B||VII-C and Appendices |D| |E]). Moreover, the discussions 
in Section IVl] allow us to design of a GRoU technique to deal with unbounded target PDFs, in Section 



IX Finally, we use the considerations and remarks in Section VI to clarify certain aspects about the 



optimality on the choice of the functions g{u). Indeed, we have deduced which is the function g{u) to 
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obtain a rectangular region Ag (see Section VI-C I. 
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Appendix 



A. Proof of the GRoU 

Given the transformation {v,u) G 



[x,z) 

V 

9{u) 



V = xg{z) 



u = z 



(91) 



z = u 

and a pair of RV's {V, U) uniformly distributed on Ag, we can write the joint PDF q{x, y) of the 
transformed RV's {X, Z) as 

1 



q{x,z) 



J-^\ for all ^<z<g-^[cp{x)\, 



(92) 



J-^ = det 



(93) 



where \Ag\ denotes the area of Ag, and J is the Jacobian of the inverse transformation, namely, 

g{z) xg{z) 
1 

Since we assume g{z) > (i.e., g increasing), then \ J^^\ = \g{z)\ = g{z) and substituting (93 1 into (92 1 
yields 



q{x,z) 



1 

0, 



g{z) for ^ < z < g '^[cp{x)], 



(94) 



otherwise. 



Hence, integrating q{x, z) w.r.t. z yields the marginal PDF of the RV X, 



]{x) 



+ 00 



l[x, z)dz 



9 ^[cpix)] 



lA 



7g{z)dz 



1 



9 ^[cp{x)] 




1 



(95) 



■5(0) 



where the first equality follows from Eq. ( [94| ) and the remaining calculations are trivial. Since we have 
also assumed ^(0) = 0, it turns out that 



q{x} 



\A, 



-p{x) = po{x). 



Hence, we have proved that a marginal PDF is exactly pq{x) oc p{x). □ 
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B. Important observation 



In the proof above, we have integrated the bidimensional pdf q{x,z) in Eq. (94 1 w.r.t. z finding the 
marginal pdf q{x) that is exactly our target pq{x). Note that the set 



An 



(v,^) € : < n < g-^ 



c p 



9{u) 



can be expressed in terms of {x,z), i.e., 

Ag = {{x,z) G . < ^ < ^-1 [c p{x)]] , 

and since g^^ is increasing and assuming p{x) decreasing (c is positive), we can rewrite it as 

1 



An 



G : < X <p"^ 



Then, if we integrate q{x,z) w.r.t. x we obtain (setting k = \/c) 



q{z) 



+00 



g(x, z)dx = 

1 

1 ^ 





g{z)x 







g{z)dx 
1 



<7(z)p-^ [k g{z)] - ^5(^)0, 



\An 



g{z)p [k g{z)] . 



(96) 



Namely, the RV Z = [/ (see Eq. (|9T])) is obtained as a transformation Z = U = g '^{\Y) = g ^(cY) 
of the RV y with pdf p^^{y) that is exactly what we anticipate in Section 



VI 



C. Other interesting observations 
It is interesting to notice that: 

• If we consider the same = n) G : < n < g~^ c p (^^^^ | (i-^-' the set Ag is defined 
in the same way) but we take x = —j^^ then we draw samples from p{—x). 

• If we consider another definition of the set = n) G : < u < g^^ c p (^~]^^ | ^i^d 
later we take x = — the set A'g is a symmetric version of Ag with respect the axis u and we 
draw samples from p{x). 

These considerations can be easily inferred from the proof above. 



D. Extension of the GRoU 

We present here a light extension of the GRoU. 
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Theorem 3: Let g{u) : [b, +00) — )• IR+ be a strictly increasing (in (6, +00)) differentiable function 
such that g{b) = and let p{x) > be a PDF known only up to a proportionality constant. Assume that 
{v,u) G is a sample drawn from the uniform distribution on the set 



■An 



{v,u) eR"^ ■.b<u< 



cp\ 



9{u) 



(97) 



where c > is a positive constant and 9 = Then x = is a sample from po{x). 

The proof is straightforward. Indeed, in this case, the development of the proof is identical yielding a 



expression similar to the Eq. ( 95 1 that becomes 



■g{z)dz 



1 



9 Mc2p(a;)] 
b 



7P{X) 



lA 



lA 



that is proportional to the target PDF (since p{x) oc Pq{x), as well). 



E. Other extension of the GRoU with a decreasing function g{u) 

We present another light extension of the GRoU where g{u) is decreasing. 

Theorem 4: Let g{u) : M+ (i.e. u < 0) be a strictly decreasing (in M"\{0} = (-00, 0)) 

differentiable function such that g{Q) = and let p{x) > be a PDF known only up to a proportionality 
constant. Assume that {v, u) € is a sample drawn from the uniform distribution on the set 



9de 



{v,u) G : g-^ 



cp\ 



9{u) 



<u<0 



(98) 



where c > is a positive constant and 5 = 3^ < 0. Then 
possibility is to define 



x 



^ is a sample from po{x). Or, another 



A' 



9de 



{v,u) £ 



cp\ 



9[u) 



< n < 



(99) 



and then take x 

It is important to note that g~^{y) : M"*" — )• M~ then g~^{y) < 0. For instance, it is possible to consider 
g{u) = /2 with u < 0, the region Ad<,c have the same form of Ag (when we use g{u) = with 

> 0) but it is symmetric the originof the axes (0, 0) , and A'g^^^ is symmetric to Ag w.r.t. the axis v. 

Finally, consider again a decreasing bounded PDF p{x) with x G M+. Then Adec can be rewritten as 
(with c = 1) 



A, 



whereas A„ can be rewritten as 

ydec 



{{v,u) : p-\g{u))g{u) < V , 



X 



u) G 



0<v< -p-\g{u))g{u)} . 



(100) 



(101) 
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Note that above —g{u) > 0. Moreover, for instance we can consider jointly Eq. ( 101 1 above and Eq. 



(51 1, and then we can write 



A = {( 



V, u) G 



■.Q<v<p-\g{u))\g{u)\] 



(102) 



where g{u) can be increasing or decreasing. Eq. (102i is clearly the expression of a transformation of 
a random variable Y with PDF p^^{y). Clearly, we can come back and obtain Eq. (99 1 or Eq. (50 1 



depending if g{u) is decreasing or increasing, respectively (note that, although we have 1(7(1*) |, however 
we also have to invert g{u); in one case the region is defined for n < and in the other for n > but 



both for V > 0). If g{u) is decreasing we obtain Eq. (99 1 and we need to set x = — whereas if g{u) 



is increasing we obtain Eq. (50 1 and we need to set x = +^^> hence finally in both we can summarize 



both cases using the following expression 



(103) 



F. Relationship between Fy^y) and Fx{x) 

The CDF -Fy(y) of RV Y ~ p^^{y) can be easily expressed as function of Fx{x) (the CDF of 
X ~ p{x)) for monotonic decreasing target pdfs pq{x) oc p{x). Indeed, we can write 

1 



K 



FxiP~Hy))+p \y)y, 



(104) 



that can be easily deduced observing Figure 18 Indeed, in Figure 18 we can see that the area -Fy(y') = 
Probjy < y'} = p~^{y)dy can be obtained as sum of the area ^ — Fx{x') = ^ — Fx{p^^{y')) and 
the rectangular area y'x' = y'p^^{y') (where we have use the relationship x' = p^^{y')). 
Observe that if we calculate the first derivative of Fyiy) in Eq. ( 104[ ), we obtain 

dFviy) 



dy 

dFriy) 
dy 

dFriy) 
dy 



dF^jp-Hy)) P-Hy) , p-]iy)y , ^-i^^^^ 



dx 



dy 



dy 



, _i. ..P ^(y) , P Hy) , -1/ N 

■p{p (y))— } — + — 1 — y + p (y)> 



dy 



dy 



(105) 



p ^{y), 



that is exactly the foreseen relationship between CDF and the corresponding PDF. 



G. Considerations about the Eq. {39) 
By definition of inverse function we have 



dig ^(cy)) = cy, 
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Fig. 18. Relationships between the CDFs of a RV Y with pdf p^^{y) and Fx{x) of X with a monotonic decreasing 

unnormalized target pdf p{x) (1/K — p{x)dx.). The area Fyiy') = Prob{y < y'} = J^^ p~^{y)dy can be obtained as 
sum of the area — Fx{x') — — Fx{p^^{y')) (depicted with oblique solid lines) and the rectangular area y'x' = y'p~^{y') 
(indicated with dashed lines). 



and then we can calculate the derivative with respect to y of both sides (using the chain rule) 

9'Hcy) 



and finally 



so that 



dg 
du 



c = c, 



du 

dg_ 

du 



dg~' 



cy 



cy 



1, 



g^^{cy) 



1 






dy 


cy 
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Acronyms and Notation 

H. Acronyms 



List of acronyms used in the text 


CDF 


Cumulative Distribution Function 


GRoU 


Generalized Ratio or Unirorms 


IID 


Independent Identically Distributed 


loD 


Inverse of Density 


MC 


Markov Chain 


MCMC 


Markov Chain Monte Carlo 


PDF 


Probability Density Fimction 


RoU 


Ratio of Uniforms 


RS 


Rejection Sampling 


RV 


Random Variable 


SMC 


Sequential Monte Carlo 


TR 


Transformed Rejection 


TRS 


Transformed Rejection Sampling 


VDR 


Vertical Density Representation 


U-GRoU 


Generalized Ratio of Uniforms for Unbounded densities 
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/. Notation 



Summary of the notation used in the text 

X Target random variable. 

Symbol used to indicate that a certain RV X is distributed according to a given proper 
normalized or unnormalized PDF. For instance, X ~ Po{x) indicates that the PDF of X is 
exactly po{x) (since p(){x) is a normalized PDF), whereas Y ~ p^^{y) indicates that the PDF 
of Y is proportional to p~^{y) (since p~^{y) is an unnormalized PDF, and the normalized 
PDF of Y is actually Kp-^{y)). 
X, X' Particular values taken by a single realization of the target RV X. 
po{x) Proper normalized PDF of the target RV X, indicated as X ~ po{x), meaning that 
Pr{X = x} = po{x). 

p{x) Proper but unnormalized PDF of the target RV X, such that pq{x) = Kp{x) for some 
K > 0. Alternatively, the relationship between pq{x) and p{x) is frequently indicated as 
Po{x) oc p{x), omitting the proportionality constant, K. The notation X ~ p{x) will also 
be used to indicate that the PDF of X is proportional to p{x), i.e. that X is distributed 
according to p{x) up to a proportionality constant if > 0, meaning that Pr{X = x} oa p{x). 
K Proportionality or normalization constant, K > Q, for the target PDF. This constant is 

independent from the value taken by the RV, x, and can be formally obtained as 





' POO 


-1 




-1 




K = 


/ p{x) dx 




/ p{x) dx 




/ P ^{y) dy 




.J —oo 











(•, ) Symbols used to indicate that an interval may be left/right open or closed. For instance, 
(a, b] indicates a right-closed interval which may be either left-closed or left-open. Similarly, 

[a, h) indicates a left-closed interval which may be either right-closed or right-open. Finally, 
(a, h) denotes an interval that may be either closed or open on both sides. 
Vx Support of the RV X, given by the range of values of X where its PDF is strictly greater 
than zero: 

Px = {a; G M : pq{x) > 0} = e M : p{x) > 0}. 



The most common supports used in the paper are Vx 
two arbitrary real numbers a, h. 



l,Vx 



and Vx = {a, h) for 
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Summary of the notation used in the text (continued) 



Fxix) Unnormalized CDF of the target RV X: 

Fx{x) = I p{x) dx. 

J —oo 

Note that we have Fx(oo) = l/K instead of Fx(oo) = 1, since we use the unnormalized 
target PDF, p{x), instead of the normalized target PDF, po{x). 
U{C) Uniform density with bounded support C, which can be a unidimensional or bidimensional 
region. In one dimension the most commonly used support is an interval starting at zero, 
[0, a], denoted by Z^([0,a]), although more general uniform distributions with other types 
of supports, such as U{px) when T>x is a bounded interval, can be considered. In two 
dimensions the most commonly used support will be the region associated to the area 
enclosed by certain PDF, e.g. U{Aq), with defined as below. It is important to remark 
that U{C) is always used to denote a proper PDF, implying that the Lebesgue measure of 
C, \C\, must be finite (i.e. that C must have a finite length for unidimensional regions or a 
finite area for bidimensional regions). 
■k{x) Proposal density from which samples can be easily drawn that is used to generate random 
variables, typically in RS and TRS algorithms. It can be normahzed or not, but always refers 
to a proper PDF. 

L Proportionality constant for the RS method, such that p{x)/tt{x) < L for some < L < oo 

and any x G Vx- 

Random variable distributed according to the inverse of the target PDF. In the sequel, this 
RV will be called the inverse target RV, and may be used to denote Y ~ PQ^{y) or (more 

often) Y ^ p^^{y)- 

y, y' Particular values taken by a single realization of the inverse target RV Y. 
PQ^{y) For monotonic target PDFs (such as exponentials or half Gaussians), proper normalized PDF 
of the inverse target RV Y, indicated as y ~ pQ^{y), meaning that Pr{y = y} = Po^iy)- 
This PDF is given by the inverse function of po{x), implying that po op^^{y) = y or 
alternatively thatpg ^opo(a;) = x, with fog{-) = f{g{-)) denoting the functional composition 
of functions / and g. 
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Summary of the notation used in the text (continued) 



Poiy) 



For monotonic target PDFs, inverse of the unnormalized PDF of the target variable, p{x), 

obtained inverting p{x) in the same way as described for Po^(y). Note that p~^{y) is a 
proper but unnormahzed PDF, since p~^{y) > for any value of y and 



dy = / p{x) dx = — , 

-OO J — OO ^ 



for K > but _ftr / 1 in general. However, unlike in the case of the target RV, the 
normalized PDF of the inverse target RV, pQ^{y), can no longer be obtained from p~^{y) 
simply multiplying Po^(y) by a normalization constant, as the support of p~^{y) can be 
different from the support of PQ^{y). Hence, a scaHng operation must be performed instead 
on p~^{y) in order to obtain pQ^{y). Formally, Y ~ PQ^iy) = P~^{y/K). Finally, we also 
note that the normalized version of p~^{y) is different from the normalized inverse target 
PDF, i.e. Kp-i(y)7^Po^(y). 

For real numbers it is used to indicate their absolute value (e.g. |a;|), whereas for a 
unidimensional or bidimensional region, C, it is used to denote its Lebesgue measure (i.e. its 
length for one-dimensional regions or its area for two-dimensional regions), e.g. |^o| = 1/^ 
is the Lebesgue measure of which is identical to the integral of p{x) over its support. 
Unnormalized generahzed inverse PDF. Given the set 

where Aq is the region enclosed by p{x), as defined below, then the unnormalized generalized 
inverse PDF is defined as 

where |.4o|y| is the Lebesgue measure of Ao\y, as defined above. 
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Summary of the notation used in the text (continued) 



Vy Support of the inverse target RV Y, given by the range of values of Y where its PDF is 
strictly greater than zero. It may be referred to F ~ Po^{y), implying that 

VY = {yeM: p^\y) > 0} = {y G M : p-^{y/K) > 0}, 

or (more often) to y ~ P~^{y), resulting in 

VY = {yeR:p-\y)>0}. 

The most common supports used in the paper are Vy = and Vy = (0, 1). Note that, 
since po{x) = Kp{x) > 0, in this case we can never have negative values of y in the support 
of Y, implying that M can never be the support of Vy. 
f{x) Invertible transformation, / : Vx — > Vz, used by the TR method to convert the target RV, 
X, unbounded and/or with infinite support, into another RV, Z, bounded and with finite 
support. 

Transformed target RV obtained using the TR method by applying the invertible 
transformation f{x) to the target RV X, i.e. Z = f{X). In some cases it is also used 
to indicate an auxiliary RV. For instance, in certain situations Z is used to denote a uniform 
RV in [0,1]. 

z, z' Particular values taken by a single realization of the transformed target RV Z. In some cases, 

it is also used to denote a single realization of ^^([0, 1]), as indicated above. 
f~^{z) Inverse of the transformation used in the TR method. 
f~^{z) Derivative of /"^{z): 



dz 



dz 



p{z) Unnormalized PDF of the transformed target RV Z = f{X), indicated as Z ^ p{z)- 
Since f{x) is an invertible function, this PDF can be expressed compactly as p{z) = 

p{r\z)) X \r\z)\. 

Vz Support of the transformed RV Z, given by the range of values of Z where its PDF is strictly 
greater than zero: 

Vz = {zeR: p{z) > 0}. 

Since the goal of the TR method is obtaining a bounded PDF with bounded support, in the 
paper we consider Vz = (0, 1) without loss of generality. 
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Summary of the notation used in the text (continued) 



Ky) 



u 



u, u 

h-Hy) 
h-Hy) 



Invertible transformation, h : Vy — )■ T>jj, used by the inverse TR method to convert the 
inverse target RV, Y, unbounded and/or with infinite support, into another RV, U, bounded 
and with finite support. 

Transformed inverse target RV obtained using the inverse TR method by applying the 
invertible transformation, h{y), to the inverse target RV, i.e. U = h{Y). In some cases, 
we may write U = U (see below) and the reason (due to the fundamental theorem of 
simulation) is given in the text. 

Particular values taken by a single realization of the transformed inverse target RV U. 
Inverse of the transformation used in the inverse TR method. 
Derivative of h~^{y): 



h-Hy) = 



dh-\y) fdh{y)\-' 



dy 



\ ^y J 



u 



q{u) Unnormalized PDF of the transformed target RV U = h{Y), indicated as Y ^ q{u). 

Since h{y) is an invertible function, this PDF can be expressed compactly as q{u) = 
p^^{h~^{u)) X \hr^{u)\ when Y ^ P~^{y), which is the case of interest considered in 
the paper 

Support of the transformed RV U, given by the range of values of U where its PDF is strictly 
greater than zero: 

Vfj = {ueR: q{u) > 0}. 

Since the goal of the TR method is obtaining a bounded PDF with bounded support, in the 
paper we consider = (0, 1) without loss of generality. 
^0 Region defined by the loD method inside which samples must be drawn uniformly in order 

to obtain samples from the target PDF, Po{x) oc p{x), and the unnormalized inverse target 
PDF, p~^{y). This region can be defined using the unnormahzed target PDF, p{x), or the 
inverse target PDF, p~^{y), as discussed in the paper. Formally, 

Ao = {{x,y) eR^ : < y < p{x)} = {iy,x) eM.^ : 0<x<p-\y)}, 

and letting {X, Y) ~ U{Ao), then X ~ p{x) and y ~ P~^iy)- 
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Summary of the notation used in the text (continued) 





Region defined by the RoU method inside which samples must be drawn uniformly in order 




to obtain samples x = v/u distributed according to the target PDF. Formally, thus region can 




be defined using the unnormalized target PDF, p{x), as 




Ar = {iv,u) e : < u < ^p{v/u)}, 




11..' / T r T T \ T 1 / A \ .1 TT" f r 1 T T /\ 

and lettmg (V, U) ~ U{Ar), then X = V/U ~ p{x). 


( \ 

9[u) 


Strictly mcreasing diiierentiable lunction on such that g[{j) = used by the GRoU 




method. Settmg g(u) = u /2 and c = 1/2 the GRoU becomes the RoU. 




Denvative of the function g\u) used by the GRoU method. 


9 [v) 


Inverse of the function g{u) used by the GRoU method. 




r^nn*;tant ^ n*;prl Hv the rrRoTT Slpttincr n(n\ — 7/^/9 anH r — 1 /9 the rrRoTT henomp*; 

V-'V./llij LclllL) O ^ Uo^U L/V Lilt/ \jAVV./\_f. O^LLlllg fJXUjJ Lt / CUIU C- A. 1 £j \XIS^ L/^C'V./lll&o 




the RoU. 


Ag 


Region defined by the GRoU method inside which samples must be drawn uniformly in 




order to obtain samples x = v/g(u) distributed according to the target PDF. Formally, this 




region can be defined using the unnormahzed target PDF, p{x), as 




Ag = {{v,u) G < u < g-\cp{v/g{u)))}, 




and lettmg (V, U) ~ U{Ag), then X = V/g{U) ~ p{x). Note that, settmg g{u) = and 




-t 1 ir\ 1.' A A 1.1 r» tt .1 ii j1 i~* x t 

c = 1/2 we obtain Ag = Ar and the GRoU method becomes the RoU. 


U, V 


T*T T 11 J 1 T* XT 1 ✓"^X^ XT jl 1 Ijl jjl * / T T T r\ ' ' f 1 l*j*1 j 1 

RVs used by the RoU and GRoU methods such that the pair [U, V) is umformly distributed 




* '1 A / C j1 X^ X T\ A / C J 1 X^ X T\ X 1 T T T T / 1\ 

mside Ar (for the RoU) or Ag (for the GRoU). In some cases, we have U = U (see above) 




1 jjl f 1 jIjI f * Ij* 

due to the fundamental theorem of simulation. 


u, u 


X^ j*l 1 j1 1 *1 l*j* f jt X^T T X T 

Particular values taken by a single realization of the RV U. 


t 

V, V 


Particular values taken by a single realization of the RV V. 


W 


Uniiorm RV used by the loD and TR methods. 


w, w' 


Particular values taken by a single realization of the uniform RV W. 




Used to denote a class function. A function /(x) is said to be of class if it is 




continuously differentiable, i.e. if f{x) is continuous, differentiable, and its derivative, f{x), 




is also a continuous function. 
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Summary of the notation used in the text (continued) 



A. 


Set of vertical asymptotes of p{x), i.e. set of points x 


G Vx where p{x) — )• oo. 


* 

X 


Used to denote the vertical asymptotes of p{x) (i.e. each of the points of X*). 


y 


Set of vertical asymptotes of p^^{y), i.e. set of points 


y e Vy where p^^{y) — )■ oo. 


* 

y 


Used to denote the vertical asymptotes of p~^{y) (i.e. 


each of the points of 3^*). 


z* 


Set of vertical asymptotes of f~^{z), i.e. set of points 


z G T>z where f~^{z) ±00. 


z* 


Used to denote the vertical asymptotes of f~^{z) (i.e. 


each of the points of Z*). 


U* 


Set of vertical asymptotes of h'~^{u), i.e. set of points 


u G Vfj where h~^{u) — >■ ±00. 


u* 


Used to denote the vertical asymptotes of h~^{u) (i.e. 


each of the points ofU*). 
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